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Hydrostatic drives as safe 
and energy saving machines 


The drive investigation method compatible 
with the diagram of power increase opposite 
to the direction of power flow 


Zygmunt Paszota, Prof. 
Gdansk University of Technology 


ABSTRACT 


Designers and manufacturers of hydrostatic drives have not a tool for precise determination 
of their energy efficiency, changing in a wide range in the driven device operating field, i.e. 
in the full range of the hydraulic motor speed and load and the working medium viscosity. 
This applies both to determining of losses and energy efficiency of displacement machines 
(pump and hydraulic motor) used in the drive system, and to losses and energy efficiency of 
the motor control system structure. A method is proposed compatible with the diagram of 
power increase opposite to the direction of power flow. That diagram replaces the Sankey 


diagram of power decrease in the direction of power flow. 


Key words: hydrostatic drive and control; evaluation of energy losses; 
method of the energy efficiency investigation; search for the energy saving solutions 


INTRODUCTION 


Hydrostatic drive and control system have been in use for 
decades in ships and ocean technology objects, mainly because 
they are safe. 

For instance, in tankers, because of the danger of explosion, 
the deck equipment — anchor windlass, automatic mooring 
winches etc. — are equipped with the hydrostatic drive. Cargo 
valves along the ship have hydrostatic remote control from the 
central console. The transport (impeller) pump system for oil 
unloading is often driven by displacement hydraulic motors 
submerged in the tanks. 

In bulk carriers the hydrostatic drive is used for large hatch 
covers, deck cranes and other deck equipment. 

Tankers and bulk carriers (and also other ships) have the 
hydrostatic drive and control systems for steering gear, bilge, 
ballast and fuel valves, cranes, hull opening covers and also 
the hydrostatic control of controllable pitch propellers in the 
ship main propulsion system. 

In tankers and bulk carriers the central hydrostatic systems 
are used of simultaneous drive and control of many devices 
with total power of several megawatts. 

A specific feature of the roll on — roll off ships are 
hydrostatically driven transport devices. A typical example is 
the stern ramp, 50 m long and 22 m wide close to the hull, and 
also internal ramps between decks. 


In the container carriers a typical example are deck cranes 
with 450kN lift capacity and 30 m or more outreach, with 
separate hydrostatic drive in three independent closed circuits 
of hoisting, turning and change of outreach. 

Fishing ships of up to 30 m in length are equipped with 
hydrostatically driven trawl, net and auxiliary winches. 
Multitask trawlers of up to 50 m in length, seiners up to 
80 m long for catching fish in shoals in the ocean depth have 
hydrostatically driven winches; the hydrostatic drive is also 
used in cranes, auxiliary winches, hatch covers and other 
centrally controlled devices. 

On the platforms, drilling, production and auxiliary ships, 
e.g. ocean tugs, the deck and technological equipment has most 
often the hydrostatic drive and control because of the danger 
of explosion. 

However, the above given examples of the use of hydrostatic 
drive in floating units do not suggest that the hydrostatically 
driven ship and ocean object equipment will be dominating 
also in the future. 

In 2006, the MacGREGOR company introduced, apart from 
deck cranes with hydrostatic drive, the first deck crane with 
electric drive [1]. In 2010 that company equipped only with 30t 
electric cranes 8 bulk carriers of 16 900 DWT capacity ordered 
in China by the Polish Polsteam owner [2]. The advantages, 
pointed out by the manufacturer, compared with hydrostatically 
driven cranes, are the following: 
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— faster and more precise positioning of the hook, reducing 
the loading time by 20% and shortening the port time, 
— the drive power smaller by 30 + 35%. 


Also proposed is the use of electric drive and control of 
hatch covers in bulk carriers as well as full electric control of 
the ro — ro ship equipment. 

The manufacturers of electrically driven ship equipment 
consider the energy efficiency of electric drive definitely higher 
than the hydrostatic drive efficiency. 

The above quoted comparison of the installed power and 
energy efficiency of electric and hydrostatic drive is made 
in the situation when designer of a hydrostatically driven 
device has no tool for precise determination of the drive 
energy efficiency changing in a wide range in the driven 
device operating field. Therefore, there is no possibility of 
precise determining of the required power and in consequence 
of a reliable comparison of the two drives. This applies to the 
losses and energy efficiency of displacement machines (pump 
and hydraulic motor) used in the drive system and also to the 
energy efficiency of the whole system as a consequence of the 
used motor speed control structure. 

In a hydrostatic drive system the mechanical energy is 
converted in the displacement pump into pressure energy of 
viscous fluid and the fluid pressure energy is converted into 
mechanical energy in the hydraulic motor. These two energy 
conversions in the displacement machines are connected with 
energy losses; the losses occur also in the conduits and in the 
motor speed throttling control assembly with a determined 
structure. 

The only way is correct assessment of the level and 
proportion of energy losses, allowing also to search for energy 
saving solutions of the pump and motor as well as the whole 
hydrostatic drive system. 


DETERMINING OF THE ENERGY LOSSES 
WITH THE ROTATIONAL HYDRAULIC 
MOTOR INVESTIGATION AS AN 
EXAMPLE 


The excerpts of the reference [19] will be quoted and 
supplemented in chapter 3 with significant findings. 

Correct evaluation of the energy behaviour of a hydraulic 
motor is an evaluation of its overall efficiency ny = f (n,,, M,,, 
v), i.e. evaluation of the overall efficiency n,, as a function 
of the motor shaft speed n,, and load M,, and of the working 
fluid viscosity v. This is also assessment of the value and 
proportions of the motor mechanical, volumetric and pressure 
losses deciding of the motor mechanical Ny» volumetric ny, 
and pressure Ny, efficiency, where the product ny = Nym Nyy ue 
determines the motor overall efficiency n,,. The energy losses 
and the corresponding efficiencies nyn» Ny, and ye should be 
determined as a function of parameters having a direct impact 
on the particular losses and efficiencies. 

However, designers and makers of rotational hydraulic 
motors and hydrostatic systems have not had so far a tool 
to determine their energy behaviour in the (0 < © M(© Me 
0<Mm<MmMna field of change of the hydraulic motor shaft 
speed and load coefficients and in the Vin < V < Vna field of 
change of the hydraulic oil viscosity. 

The rotational motor producers make erroneous routine 
evaluations of the following energy efficiencies and work 
parameters: 

— the motor overall efficiency ną as a function of the shaft 

speed n,, and motor pressure decrease Ap,, (e.g. [11 + 

14), 
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— the motor overall efficiency ną as a product of the volumetric 
efficiency n,,, and the so called motor „mechanical 
— hydraulic efficiency” nyp all the three determined as 
a function of the same parameters (e.g. [11, 12]), 

— motor shaft speed n, as a function of the motor capacity Q,, 
and volumetric efficiency n,,., determined in turn as a function 
of the motor pressure decrease Ap,, (e.g. [11, 13]), 

— motor shaft torque M,, as a function of pressure decrease 
Ap,, and the so called „mechanical — hydraulic efficiency” 
Nsimn OF the motor (e.g. [11 + 14]), 

— motor capacity Q„ as a function of the shaft speed n, and 
volumetric efficiency Ny, determined in turn as a function 
of the motor pressure decrease Ap,, (e.g. [11 + 13]), 

— motor shaft useful power P,,, as a function of the motor 
capacity Q,, and pressure decrease Ap,, and as a function 
of the motor overall efficiency ny determined in turn as 
a function of the motor shaft speed n, and pressure decrease 


Ap,, (e.g. [11 + 13). 


The hydraulic motor researchers evaluate in a wrong way 
the losses arising in the motor: 

— the motor torque M „m of mechanical losses as a function 
of the motor pressure decrease Ap,, and shaft speed n,» 

— asum of torque M,,, of mechanical losses and the so called 
„torque of pressure losses” (resulting from the pressure 
losses Ap,, in the motor) — as a function of motor pressure 
decrease Ap,, and the shaft speed n» 

— the intensity Q„, of volumetric losses in the motor as 
a function of the motor pressure decrease Ap,, (or as 
a function of the motor shaft torque M,,) and as a function 
of the motor shaft speed ny- 


The evaluation methods of the energy losses and 
efficiency of the rotational hydraulic motors, used so far in 
the scientific research and in the industrial practice, give 
wrong results because the parameters that the losses and 
efficiencies are a function of are themselves dependent of 
those losses. 

There is very little information of the motor makers 
presenting properly the motor overall efficiency ny = f (n,,, 
M,,) as a function of the motor shaft speed n,, and torque M,, 
at a specified fluid viscosity v and presenting the impact of 
viscosity v on the overall efficiency n,, (e.g. [10]). 

It is a common deficiency that no information is given 
about the dependence of the motor mechanical, volumetric 
and pressure losses on the kinematic viscosity v of the working 
fluid used in the hydrostatic drive system. 

The fundamental reason of the erroneous evaluations are 
commonly accepted views on the research methodology and 
on the method of determining the energy losses in pumps and 
in hydraulic motors. That method is based, among others, on 
the traditional reading of the energy balance of a hydrostatic 
drive system from the Sankey diagram [3 + 9]. The present 
unsatisfactory state is also an effect of using simplified 
evaluations of the relations of particular losses to the motor or 
pump working parameters and to the working fluid viscosity. 

In reference [19] the work parameters, power, losses and 
energy efficiencies in a rotational hydraulic motor are defined 
and described by means of 46 expressions and also their 
complex interrelations are shown. 

The end formula (46) presents the motor efficiency n,, as 
a product of efficiencies described by means of losses and 
parameters deciding about those efficiencies, and the pressure 
decrease Ap,,, indicated in the motor working chambers and 
capacity Q„ in the motor chambers are also functions of the 
losses: 


R= Py Oy; 


direction 
of power flow 


viscosity » of the working fluid 


inlet channel (py ) 


u= far Rae = (Pur Puz) UFAN AE Ris + Ay 


"shaft — working 
chambers” assembly 


= 


direction 
of power flow 


shaft speed ny (c,) 
shafi torque My 


Q70! 
direction L | 
lof power flow 


outlet channel (Apyo2) / 


Ri Puu; outlet pressure py.*0 


P =R u + Aim t+ Afu +tAfip 


direction 
of power flow 


"shaft — working 
chambers assembly 


m 


direction 
of power flow 


outlet channel, 


direction 
of power flow 


Fig. 1. Diagram of power increase in a hydraulic motor opposite to the direction of power flow, 
replacing the Sankey diagram of power decrease in the direction of power flow 


Nm = f(My.0y-V) = NmmNm p= 


__ Mw qm m AP wi 
My +Myum Ant Am + Quy APmi + AP op 


where: 
Mum = f( My, 2w V) [equation (2)], 


Qus =f (APup By V) [equation (14)], 


(46) 


AP vip =f( Qis v) [equation (25)], 
with: 
APyi = 211M y + Myn) [equation (6)] 
qm 
and 


Qui = dmr Bu T Qu, [equation (16)]. 


Decrease Ap „; of pressure indicated in the working chambers 
(equation (6)) is a function of the loading torque M,, and torque 
M um of mechanical losses in the „shaft — working chambers” 
assembly. Capacity Q,, in the motor channels (equation (16)) 
is a function of the motor shaft speed n,, and the intensity Qu, 


of volumetric losses in the working chambers. 


POWER IN THE DRIVE SYSTEM, LOSSES 
AND ENERGY EFFICIENCY, SYSTEM 
OPERATING FIELD 


Diagram of power increase in a power transmission 
system 
opposite to the direction of power flow 
replaces 
the Sankey diagram of power decrease 
in the direction of power flow. 


The Sankey diagram is the fundamental cause of incorrect 
evaluation of losses in a drive system. 

The Sankey diagram connected with a drive system suggests 
an evaluation defining the output power of the system as 


a difference between the system input power and sum of the 
powers of losses occurring in the system. Therefore the energy 
losses in the system are evaluated as a function of parameters 
describing the input power. 

In accordance with the proposed diagram of power increase 
in a power transmission system opposite to the direction of 
power flow, the system input power is a sum of the system 
output power and powers of losses occurring in the system. 
The energy losses in the system are evaluated as a function 
of parameters describing the system output power. 

Fig. | illustrates this principle by an example of power flow 
in a hydraulic motor. 

The power P „, consumed by the hydraulic motor is a sum 
of motor shaft useful power P,,, and powers of three different 
energy losses in the motor. The losses occur in series increasing 
power in the opposite direction to the direction of power flow. 
In effect, the power in the motor increases from the shaft 
useful power P „ to the working fluid power P,,, consumed 
by the motor: 


Pi = Piia + AP im + AP w + AP wp 


Mechanical losses (and power AP „m of mechanical losses) 
occur in the „shaft — working chambers” assembly. 

Volumetric losses (and power AP „ of volumetric losses) 
occur in the working chambers. 

Pressure losses (and power AP,,, of pressure losses) occur 
in the motor channels. 

Formula (46) (equation (2)) shows a direct dependence of 
the torque M,,,, of mechanical losses in the „shaft — working 
chambers” assembly on the torque M,, and on the motor shaft 
rotational speed n,, as well as on the working fluid viscosity v. 

Formula (46) (equations (14, 6)) presents a complex 
dependence of the intensity Qu, of volumetric losses in the 
working chambers on the shaft loading torque M,, and on 
the torque M of mechanical losses in the „shaft — working 
chambers” assembly (decrease Ap,,, of pressure indicated in 
the working chambers depends on M,, and M and has direct 
impact on Q,,,) and also on the shaft speed n,, (influencing in 
a diversified way the torque M „„ of mechanical losses (2) and 
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intensity Q,,, of volumetric losses (14)). The intensity Qw, 
of volumetric losses depends on a diversified impact of the 
working fluid viscosity v: indirectly by impact of v on the torque 
M,, of mechanical losses in the „shaft — working chambers” 
assembly (2) and directly by impact of v on the intensity Qu, 
of losses in the working chambers (14). 

Formula (46) (equations (25, 16, 14, 6)) presents also 
a complex dependence of losses Ap„, of pressure in the motor 
channels on the motor shaft rotational speed n, and on intensity 
Qu, of volumetric losses in the motor working chambers. 
The intensity Q,,, of volumetric losses influences the motor 
capacity Q„ (equation (16)) and at the same time Q,,, depends 
in a complex way on the shaft loading torque M,, and on the 
torque M,,,, of mechanical losses in the motor „shaft — working 
chambers” assembly. Pressure losses Ap, , in the motor channels 
are also dependent on the diversified impact of the working 
fluid viscosity v: indirectly by impact of v on the torque M,,,_ of 
mechanical losses in the „shaft — working chambers” assembly 
(2) and by impact of v on the intensity Q„, of volumetric losses 
in the working chambers (14) and directly by impact of v on 
the losses Ap,,, of pressure in the channels (25). 

Contrary to the commonly used, both by manufacturers 
and researchers, methods of evaluation of the hydraulic motor 
losses, it is unacceptable to create a ,,sum” of the torque 
M um of mechanical losses in the „shaft — working chambers” 
assembly and the ,,torque” of pressure losses Ap,,_ in the motor 
channels, and also such a ,,sum” must not be evaluated as 
directly dependent on the same chosen parameters, because 
those losses are of different character and depend on different 
parameters: 


Mim = f (M,e Dy V) 
Apu, = £ (Quo V) 


The impact of hydraulic oil viscosity v on the energy 
losses in a hydraulic motor, i.e. on: 
— torque M,,,, of mechanical losses in a „shaft — working 
chambers” assembly, 
— intensity Q„, of volumetric losses in the working 
chambers, 
— pressure losses Ap,,, in the conduits, 
is diversified. 


Dependence of the motor losses on the hydraulic oil 
viscosity v should be presented in expressions describing also 
the dependence of those losses on other parameters which 
influence them directly: 


Mim > f (M,e Dap V) 
Quy = f (AP is Aw V) 
Apu, = £ (Qo V) 


The motor overall efficiency n,,, as a function of M,,, n,, 
and v, is a product of nyn» Nuy and n,,, efficiencies: 


m 7 f Mye Du Y) = 


Pw | Mm@m _ 2H Myn N NN 
= = = “~~ IMm 'Mv 'Mp 
Pue = =APmQu Ap w@Qmu 
where: 
Pu, ~ the motor useful power, 
Pi. -~ the motor consumed power. 


Each of the three efficiencies as a factor in the product 
describing the overall efficiency, is evaluated as a function 
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of parameters directly influencing the respective losses and 
a function of parameter to which the losses are ,,added”: 
— motor mechanical efficiency nyn: 


n — Pin = Myu pe 
Mm ` B 7 
Pu Mu +Mum)@m 
2IIM, n Mu 

=—_—__—___ MM f Muny V) 

2M + My en ay Mu +M iini 

where: 

Pwa — the power indicated in the motor working 

chambers: 


Pui = APmiAm Dm= 
= (Mu +Mum) Ou = 2 (Mu + Main OM 


— motor volumetric efficiency n,, : 


Mic Pui _ AP MidaieOm — APmiAmAM _ 
My = = = = 
Pui  APmi AmNm Qum)  APmiQu 

n n 
2 Jmm Ame =f (AP ny V) 
dm Dm + Qw Qu 
where 
Pie the power consumed in the motor working 
chambers: 


Puci = AP i i'm + Qu) = APmiQui 


— motor pressure efficiency nyp: 


n Pmi AP i Amm + Quy) _ AP MiQu _ 
= Pue (AP vi + AP Mp MAmAmt Qm) APMQm 
AP wi AP wi 
= = = f (APm Que Y) 
AP wi + AP ip APu =r 


In order to present the motor volumetric efficiency iig 
a factor in the nym Nu, Nip product describing The i.e. to al 
Nyy as a complex dependence on the (M,,, ny» V) parametrs 
describing n, and dependent on the ie Bee losses, the 
intensity Qu = f (Ap,,, Dw V) of volumetric losses in the 
working chambers should be determined with: 


_ 2i (My T Mum) 
Aut 


of mechanical losses in the ,,shaft 
an f (Myu Dy V) 


Ap Mi — 


and with torque M ym 
— working chambers” assembly as an M,, 
function. 

In order to present the motor pressure efficiency Np 48 
a factor in the nym Ny, un product pip ra i.e. to present 
Tmp 28 a complex dependence on the (M,,, OV) parameters 
describing Ny and dependent on the ene ‘and volumetric 
losses in the t motor, the pressure losses Ap,,, = f (Q,,, V) in the 
channels must be determined with: 


Qu = dm By t Quv 
then intensity Qu, = f (Ap,,,, Dyp V) of volumteric losses in 
the working chambers must be determined with: 


ap, = My + Mun) 
Ma qm 
and the torque M „„ of mechanical losses in the „shaft — working 


chambers” assembly must be determined as an M,, 
np V) function. 


m =f Mw 


The characteristic of the hydraulic motor overall 
efficiency ny = f (M,, n,, Y) presents a complex picture 
as a product Nym Nyy Nm, Of three efficiencies correctly 
described by: 

mechanical efficiency Nym = f (Mpp Dw V) 

— volumetric efficiency n,,,=f (Ap, My) V) 
— and pressure efficiency ny = f (APup Qu V)- 


The w of the hydraulic motor overall efficiency 
Nyu =f (M,e Ny V) must be supplemented by assessment of the 
motor ean field (0 < OMO Mma s 0<Mm(Mmna in the 
hydrostatic drive system, i.e. assessment of the range of ny 
and M,, (®,, and My parameters. 

In the motor (and the drive system) (0<O\(Oy,,. 
0<My Mm, ) operating field its current speed n, (0, *) 
and load M My are an effect of the demand of the motor 
(system) driven machine (device) and are independent of 
the losses in the hydraulic motor and in the motor driving 
hydrostatic system. 

However, the hydraulic motor (system) operating field 
limit values Nymay (O Mma) and Muma (Mmaa) depend on the 
maximum capacity of the motor driving system. The values 

Aymar (2Mmax) 2nd Mimar (Munad determine simultaneously 
the corresponding motor energy efficiency n,, and the overall 
system efficiency n. 

The limit parameters of the hydrostatic drive system 
operation result from the pump theoretical capacity Q,, and 
the system nominal working pressure p, as well as from the 
energy losses in the hydraulic motor, pump and conduits and 
also losses in the motor speed throttling control assembly (if 
is installed). Therefore, the limit values n,,, (Minox) and 

Miis (Matera are also dependent on the working fluid viscosity 
v changing 1 in the Vmin = V Š Ving, range. The hydrostatically 
driven hydraulic motor operating field is also influenced by the 
n,= f(M,) characteristic of the electric or combustion motor in 
the pump driving system. 

Evaluation of the hydraulic motor overall efficiency 
Ny = f (M,, ny, V) as a product Nym Nys Nis of three 
motor efficiencies should be performed by means of the 
mathematical models of losses and efficiencies, where the 
defined coefficients k, of energy losses in the motor and in 
the motor drive system are used. 

Evaluation of the hydraulic motor energy efficiency 
is performed together with evaluation of the efficiency 
of a hydrostatic drive system where the motor is used 
[including also the energy efficiency of pump, conduits and 
the motor speed throttling control assembly (if it is used)]. 

In the proposed metod, based on the mathematical models of 
losses, each kind of energy losses is a function of parameters 
directly influencing the losses and independent of those 
losses. 

Evaluated are the values of the k, coefficients of losses, 
relating the mechanical, volumetric and pressure losses in 
the hydraulic motor, pump and other system elements to 
the reference values of driving system nominal pressure p, 
theoretical capacity Q, of the system driving pump, theoretical 
speed n, ,, and theoretical torque M,,, of the motor shaft as well as 
the theoretical speed n,, and theoretical torque M,, of the pump 
shaft. The k, coefficients are determined at the hydraulic oil 
reference viscosity v,. At the same time the impact is determined 
of the viscosity ratio v/v, in the v in <v < Vna range on the k, 
coefficient values. 

The method allows to evaluate the values and proportions 
of mechanical, volumetric and pressure losses in the motor, 
pump, conduits and in the throttling assembly (if installed), as 
well as the dependence on the hydraulic oil viscosity v. 


The energy investigations of a pump and hydraulic 
motor as independent machines are limited to determination 
of the k, coefficients of losses in them. The energy efficiency 
characteristics of those machines are determined in parallel 
with efficiency evaluation of the system where they are 
used. 

The knowledge of k, coefficients of the mechanical, 
volumetric and pressure losses in the drive system elements 
allows to obtain, with the numerical method, the characteristics 
of the hydrostatic drive system overall efficiency, pump 
efficiency, hydraulic motor efficiency, conduit efficiency and 
the motor speed throttling control assembly (if it is used) 
efficiency in the motor (system) (0 <p Mnp OS Mud Mma 
operating field at a selected ratio v/v, of the hydraulic oil 
viscosity and the reference viscosity. 

Characteristics of energy efficiency of the pump and 
hydraulic motor with the determined constant coefficients k, 
of losses and also of efficiency of the conduits are different 
in hydrostatic systems with different motor speed control 
structures. 

Characteristics of the overall efficiency of elements used 
in a hydrostatic drive system: pump n, hydraulic motor 
Ny conduits n, and throttling control assembly (if it is 
used) structural efficiency n, are defined as functions of 
the hydraulic motor (system) “speed coefficient © „and load 
coefficient M,, and the hydraulic oil viscosity ratio V/V 

At the same time the hydraulic motor (0 < © OMO Mo 
0< Mu (Muna ) operating field in the hydrostatic drive 
system is determined at the selected ratio v/v, of the 
hydraulic oil viscosity to the reference viscosity. 

The method is precise by definition and simple in use. 
It simplifies the laboratory investigation of pumps, hydraulic 
motors and hydrostatic drive systems. It allows to seek for 
energy saving solutions of pumps and hydraulic motors. It 
allows also to evaluate the overall energy efficiency of the drive 
and finding energy saving hydrostatic drive system structures. 


CONCLUSIONS 


1. The methods of investigating the rotational hydraulic motor 
losses and energy efficiency, used in scientific research and 
industrial practice, give incorrect evaluations, because: 

— the losses and efficiencies are evaluated as functions of 
parameters which depend on those losses or as functions 
of parameters, which have no direct impact on the 
losses, 

— the mechanical, volumetric and pressure losses and 
the corresponding energy efficiencies are presented as 
dependent on the same parameters, although in fact each 
of those losses is a function of different parameters and 
is a different function of the hydraulic oil viscosity v. 

2. The basic reason for the incorrect evaluation are common 
views on the methods of carrying out the investigations and 
on the method of defining the energy losses in pumps and 
hydraulic motors. The method is based on the traditional 
understanding of the hydrostatic drive system energy 
balance illustrated by the Sankey diagram of the power 
decrease in the system in the direction of power flow 
[3 + 6, 15, 17, 19, 21, 22]. 

3. Diagram of power increase in a power transmission 
system opposite to the direction of power flow replaces 
the Sankey diagram of power decrease in the direction 
of power flow. 

4.The characteristic of the hydraulic motor overall efficiency 
Ny = £(M,, Ny, V) as a product nym Nyy, Ny, Of three efficiencies 
correctly described by the characteristics of the mechanical 
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10. 


11. 


12. 


8 


efficiency Num = f (M,, ny, V), volumetric efficiency 
7 = f( APs Dy V) and | se efficiency nw = f ( Ap,» 
ww V), is very complex. 
Se of ag hydraulic motor overall efficiency 
Ny = f (My, 2y V) must be supplemented by evaluation 
of the motor ti <Dy < OM inex ,0<SMy < MMag) Operating 
field in the hydrostatic drive! system, i.e. of the ny and M,, 
(© p and M,,) range. The hydraulic motor operating field 
limit values nyimas O Mma) and Muma M (Mma) depend on 
the maximum capacity ‘of the motor driving hydrostatic 
system. That capacity is a result of the pump theoretical 
capacity Q,,, system nominal pressure p and energy losses 
in all the system elements, which depend also on the used 
motor speed control structure. 
ates of the hydraulic motor overall efficiency 
= f (M,, ny, V) as a product nym Nm Nip of the three 
me efficiency should be based on the mathematical 
models of losses and efficiencies where coefficients k, of 
the energy losses in the motor and in the motor driving 
system are used. Proposed evaluation of the hydraulic 
motor efficiency is performed together with evaluation of 
the energy efficiency of pump, efficiency of conduits and 
efficiency of the throttling control assembly (if it is used) 
as well as the overall efficiency of the hydrostatic drive 
system, where the pump, motor, conduits and control system 
assembly are applied. 
Evaluated are the values of k, coefficients of losses 
relating the mechanical, volumetric and pressure losses 
in the hydraulic motor, pump and the whole system to the 
reference values: theoretical capacity Q, of the system 
driving pump, system nominal pressure p,, motor shaft 
theoretical speed n,,, and theoretical torque M,,, as well as 
the theoretical speed n,, and theoretical torque M,,, of the 
pump shaft. The k, coefficients of losses are determined at 
the hydraulic oil reference viscosity. At the same time the 
impact of the oil viscosity ratio v/v, (viscosity in the v_ 
<v<v,,,, range) impact on the values of k, coefficients of 
losses is determined. 
The method allows to evaluate the values and proportions 
of the mechanical, volumetric and pressure losses in the 
hydraulic motor, pump and the system and also their 
dependence on the hydraulic oil viscosity v. 
Energy investigations of the pump and hydraulic motor as 
independent machines are limited to determination of the k, 
coefficients of losses. The energy efficiency characteristics 
of those machines are determined in parallel with evaluation 
of efficiency of a hydrostatic system where the machines 
are used. 
Knowledge of the values of k, coefficients of losses 
allows to obtain, by a numerical metod, of the hydraulic 
motor efficiency ny = f(m» My, WV,), pump a 
Np= f(O My, WV,), conduit efficiency 1), =f(Oy, Mys WV,), 
throttling p assembly (if used) structural efficiency n; 
=f(©m Mu Wa)and = hydrostatic drive system overall 
efficiency n = f@y, M m» W/V,) in the hydraulic motor 
(system) (0 < Om < Mna 0<My, < MM. ) operating field 
and at the selected ratio vy, of hydraulic ‘oil viscosity. 
The method is precise by definition and simple in use. It 
simplifies the laboratory investigations of pumps, hydraulic 
motors and hydrostatic drive systems. It allows to evaluate 
the characteristics of the system overall energy efficiency 
as well as to seek for the energy saving hydrostatic drive 
structures. 
The method is a proposition for the pump, hydraulic 
motor, hydrostatic drive and control system designers and 
manufacturers. It is in the process of development. 
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Three-dimensional cfd analysis to study 
the thrust and efficiency of a biologically- 
inspired marine propulsor 


Lamas M. I., Assoc. Prof. 
Rodriguez J. D., Assoc. Prof. 
Rodriguez C. G., M. Sc. 
Gonzalez P. B., Prof. 

Escola Universitaria Politécnica. 
Universidade da Corufia 


ABSTRACT 


Aquatic animals, which are the result of many millions of years of evolutionary optimization, are very 
quick, efficient, robust, and versatile. Accordingly, biologically-inspired mechanisms which emulate the 
movement of animals have recently become very popular. For the efficient design of a propulsion system 
it is very important to analyze the fluid flow in detail. CFD (Computational Fluid Dynamics) has become 
a powerful technique to understand the phenomena because it gives extensive information about the fluid 
flow characteristics. 
In the present work, a propulsion system consisting of an undulating fin which emulates the fish swimming 
was built. In order to optimize the mechanism, several undulating configurations were studied using 
a 3D turbulent CFD model. The thrust, drag, efficiency and hydrodynamic characteristics were analyzed. 
Furthermore, it was shown that the efficiency and thrust depend strongly on the oscillation frequency, 
amplitude and wavelength. 
In order to validate this CFD model, the numerically obtained thrust was successfully compared with 
experimental results from the laboratory mechanism. 


Keywords: fish swimming; marine propulsion; undulating fin; biomimetic; computational fluid dynamics; CFD 


INTRODUCTION 


Biomimetic is an emerging field which employs the 
principles of living organisms to derive man-made mechanisms 
which are capable of emulating the efficient movement 
of animals. In the field of marine propulsion, undulating 
fish-like mechanisms offers many advantages over rotary 
propellers. What makes the high efficiency and high thrust of 
the biologically movement is the manner in which it works. 
Behind the fin, several vortices are formed with a rotation 
direction compatible with the direction of thrust, producing 
an efficient jet. On the contrary, a propeller generates a long 
jet that rotates in the direction perpendicular to the motion, 
wasting a considerable amount of power that goes into rotating 
this jet. 

Under the same conditions, a fish consumes much less 
energy to displace than a rigid body because the motion of the 
fish contributes to a reduced drag force and increased propulsive 
efficiency by reducing separation and suppressing turbulence 
(Shen et al., 2003). Recent works have shown that undulating 
fish movement is very efficient compared to rigid bodies (Barret 
et al., 1999, Liu and Kawachi, 1999, Triantafyllou, 1995) and 
there are several researchers interested in fish-like underwater 
and surface ships. For example, Barret et al. (1999) constructed 
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a fish-like mechanism and studied it on a laboratory in order 
to compare the undulating fish movement with a rigid wall 
movement. They concluded that the power required to propel 
an actively swimming fish-like body is significantly smaller 
than the power needed to tow the body straight and rigid and 
they obtained drag reductions up to 70% for the undulating 
fish movement. Zhang ef al. (2008) investigated the bionic 
neutral network control method for fish-robots and developed 
a 2D laminar CFD model to study several amplitude envelopes. 
Young-hua et al. (2007) designed an environment-friendly 
propulsion system mimicking undulating fins and developed 
a two-dimensional CFD method to study the unsteady flow 
around the fin. They focused its study on studying the amplitude 
configuration. Bozkurtass et al. (2008) constructed a propulsor 
for an autonomous underwater vehicle based on the mechanical 
design and performance of a sunfish pectoral fin, and developed 
a three-dimensional numerical model. Low (2008) constructed 
an underwater vehicle using a fin-like mechanism based on 
a series of connecting linkages which produce undulations 
similar to those produced by the fin rays. Clark and Smith 
(2007) built a flexible fin which replicates some features of 
the pectoral fin of a batoid fish. 

Although these studies have revealed many key points of 
undulating propulsion, much is still unknown regarding the 


behavior of the unsteady fluid flow. Recently, CFD has become 
a powerful tool to study these phenomena. By solving the 
Navier-Stokes equations, CFD provides information on the 
pressure and velocity fields. It is very useful to complement 
traditional testing and experimentation, providing faster time 
and lower costs. In this work, an undulating propulsor which 
emulates the high performance and maneuvering of fish 
swimming was built. In order to analyze the flow around the 
fin, a tri-dimensional turbulent CFD model was developed. 
A previously 2D CFD model was already developed elsewhere 
(Lamas et al., 2010). In the present paper, the 2D model was 
extended to the more realistic 3D situation. This numerical 
model was satisfactory compared with the experimental data 
obtained from the experimental prototype. Once validated, the 
CFD analysis was employed to predict how the performance 
varies with different parameters. 

This paper begins with a introduction of the experimental 
prototype and the kinematics of locomotion. The numerical 
method and results are then briefly described and contrasted 
with experimental data. Finally, the conclusions of this work 
and the areas for future research are presented. 


PROBLEM DEFINITION 


Design and kinematics 


The experimental ship is shown in Fig. 1.a. The propulsive 
fin, which is shown in Fig. 1.b, consists of a flexible surface of 
0.52 m length and 0.2 m width. 


b) 
y 
Flow direction 
— — 
——> —_ 
ed —S>" 


L 


Fig. 1. a) experimental mechanism, b) detail of the fin 


Fig. 2 provides the mechanism of the undulating fin, 
consisting of eight rigid segments, each capable of relative 
rotation with respect to its neighboring links. The movement 
of each segment is produced by a crankshaft situated above 
the fin. The power is supplied by an electric motor and the 
oscillation frequency is controlled by a variable-frequency 
drive control. 

Every module is attached to a connecting rod on both 
sides. The position of every segment is defined by the rotation 
of the crankshaft, which moves the rods, the connecting rods, 
and finally the segment. In order to investigate what kind of 


Fig. 2. Detail of the undulating surface and mechanical transmission 
of movement 


amplitude envelope provides the maximum efficiency, the 
mechanism was designed to accommodate several centimeters 
peak-to-peak amplitude. More details of the system were given 
elsewhere (Rodriguez, 2009). 

The high thrust and efficiency are not the only advantages of 
this system. It is also reversible, i.e., it has the same efficiency 
either operating forward or backward. This makes it ideal for 
vehicles that require high maneuverability. Other advantage is 
that the use of several properly positioned fins can minimize 
vibrations. 


Governing equations 


The fluid flow around this mechanism was studied 
numerically. In order to simplify the model, the flow was 
considered incompressible and isothermal. The unsteady 
turbulent flow was simulated by the RANS (Reynolds-averaged 
Navier-Stokes) equations. In Cartesian tensor form, the mass 
conservation equation is given by: 

õp oO 
Bt Ws (pu;) =0 (1) 


and the momentum conservation equation: 
£ (pu) z (Puu)= 

= Op + ann + 0 (puu) 

ôx a, d 


Since the fluid was treated as Newtonian, the stress tensor 
components are given by: 
ou, ĉu; 2. ðu 
a aa 6) 
j OX, Xk 
A common method to model the Reynolds stresses, — pu; u; 


is the Boussinesq hypothesis to relate the Reynolds stresses to 
the mean velocity gradients: 


(2) 


—ouu.=u,/—+—|]-Zlok+u, — |8. (4 
pu; ne a f Hy | 1J ( ) 
where ò is the Kronecker delta (6,,= 1 if i =j and 6,=0 if 


i # J) which is included to make the formula applicable to 
the normal Reynolds stresses for which i = j (Versteeg and 
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Malalasekera, 2007) and 1, is the turbulent viscosity. Since the 
turbulent viscosity was modeled by the k-e model, two transport 
equations are added, one for the turbulence kinetic energy k (eq. 
5) anda further one for its dissipation rate e (eq. 6) respectively. 
Unfortunately, no single turbulence model is universally 
accepted. The reasons why the k-e model was applied in this 
work are that it is robust, computationally economical and 
reasonably E for a wide o of turbulent flows. 


5 (PK) + (pku; )= 


mn (Peu) = a ut) 3 + 
Ox; e j (6) 
€ 
C k = (G, + Gz; G,)- Ca P PE 
In the above equations, G, represents the generation of 
turbulence kinetic energy due to the mean velocity gradients; 
G, is the generation of turbulence kinetic energy due to 
buoyancy; Y,, represents the contribution of the fluctuating 
dilatation in compressible turbulence to the overall dissipation 
rate. C Cie C,,, Cso 6, and o, are constants with the standard 
k-e default values C, = 0.09; C,, = 1.44; C,,= 1.92; o,= 1.0 
and o,= 1.3. Details are given in (Fluent Inc., 2006). 
The turbulent viscosity, u, Is computed by combining k 
and e as follows: 
. 
=pC,— (7) 


NUMERICAL IMPLEMENTATION 
Computational mesh 


In order to implement the movement of the fin, a moving 
mesh which regenerates with time was employed. The 
computational domain is shown in Fig. 3, in which the front 
surface and the internal mesh were omitted for clarity. It was 
discretized with 600000 tetrahedral elements. The size mesh 
was refined near the fin and in the wake region to resolve the 
shed vortices as they travel downstream. Due to the symmetry 
of the problem, only half of the geometry was modeled. 


Fig. 3. Computational domain 
Boundary Conditions 
The downstream boundary (right surface) was modeled 


as a pressure outlet and the upstream boundary condition 
(left surface) was modeled as an uniform velocity inlet 
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corresponding to the free-stream velocity. The fin surface was 
treated as a no-slip moving wall. The top and bottom sides were 
modeled as no slip walls and the front surface was treated as 
a symmetry plane. 


Numerical Methodology 


The numerical solution was performed using the commercial 
software Ansys Fluent 6.3, which is based on the finite volume 
method. Pressure-velocity coupling was achieved using 
the PISO algorithm. A second order scheme was used for 
discretization of the pressure and momentum equations. The 
time integration was achieved using a first order fully implicit 
scheme with a constant time step of 0.001 s. 

Extensive convergence checks were taken in order to 
confirm that the computed results are independent of the time 
step, the grid size and the computational domain size. 


CALCULATION OF FORCES 
AND EFFICIENCY 


As the fin moves through the water, a force is produced 
along the x, y and z directions. Particularly, a pressure force 
is caused due to pressure difference and a friction force is 
caused due to water viscidity. The components of the forces 
acting on the body, F, F, and F, can be evaluated by integrating 
the projection of the pressure ‘and the shear stress in the x, y 
and z directions respectively. The total component along the 
x direction was computed by adding the component x of the 
pressure and viscous forces acting on the body as follows: 


F= P +F. (8) 


where: 
F — the pressure force 


px 


F  — the viscous force in the x direction. 
The pressure force along the x axis is given by: 
F,=— [pnas (9) 
A 


where: 
n, — the x component of the unit normal vector on dS. 


The viscous force along the x axis is given by: 


F= frgndS (10) 
A 
where: 
t. — the viscous stress tensor. 


The total power required for the propulsive motion of 
the fin (P,) consists of two components, P ,= P, + P, One 
is the swimming power (P) required to produce the vertical 
oscillations of the fin and the other is the useful power (P,,), 
required to thrust. The swimming power was computed as: 


P= fp AY, at ds (11) 
dt 

where p is the pressure of the surface of the plate and dy „dt 

refers to the y-component of the velocity at the surface of 

the fish. For a sinusoidal wavy the motion of the fin and the 

y-component of the velocity are given by Ec. (12) and (13) 


respectively. 
: t 
y = A(x): sin [x = — + (12) 


Hoa. Us - cos] 27 cea (13) 
dt T A T 
where: 
A(x) — represents the amplitude, 
À — the wavelength, 
T  — the period. 


The useful power was computed as the time-average thrust 
multiplied by the free-stream velocity: 


PeP U (14) 


where the minus sign indicates that the velocity of swimming 
is reverse to that of the stream. 

The efficiency was computed as the ratio between the useful 
power and the total power. Obviously, the efficiency is always 
less than one because of losses due to wayward vortices and 
turbulence. 


RESULTS AND DISCUSSION 
Configurations Studied 


The motion of the fin is in the form of a streamwise 
travelling wave. As the fluid moves along the surface, the 
undulations from the anterior to the posterior of the fin produce 
thrust. The amplitude may vary along the length, and its design 
is very important. For this reason, several CFD analysis were 
carried out to assess the most suitable configuration. Three cases 
were studied: amplitude increasing from the anterior part to the 
posterior (Fig. 4a), constant amplitude (Fig. 4b) and amplitude 
increasing from the anterior part to the mid part and decreasing 
towards the posterior (Fig. 4c). 

Ay en aS 


Fig. 4. Amplitude configurations. a) amplitude increasing from the anterior 
part to the posterior, b) constant amplitude, c) amplitude increasing from 
the anterior part to the mid part and decreasing towards the posterior 


Hydrodynamic forces and efficiency 


Fig. 5a represents the time history of the instantaneous 
pressure force, viscous force and average total force along the 
cycle at U „= 0; f= 10 Hz; A „= 0.02 mand A= 1 m. The zero 
free-stream velocity corresponds to a experimental situation 
in which the ship is tied to a fixed point by means of a rope, 
Fig. 5b. The configuration studied in Fig. 5a corresponds to (a), 
i.e., amplitude increasing from the anterior part to the posterior, 
and the other two amplitude configurations presented a similar 
tendency, so they are not plotted again. 


In order to reach a situation which periodically repeats, i.e., 
a quasy-steady state, it was necessary to study a long enough 
interval of time. For the cases studied it was verified that this 
state is achieved after approximately twenty fin oscillation 
periods. For this reason, all the results carried out in the present 
paper correspond to the 20* period of time. 
a) 20 


o-—- Pressure force 
Viscous force 
Average total force 


15 


10 


Force [N] 


0 0.2 0.4 0.6 0.8 1 
Time (t/T) [-] 


Fig. 5. a) time history of the instantaneous thrust and drag and 
average total force at U, = 0; f= 10 Hz; A „„ = 0.02 m and) = 1 m. 
Configuration (a), amplitude increasing from the anterior part to the 

posterior; b) equivalent experimental situation 
(pulling against a fixed point) 


The forces shown in Fig. 5a are the ones that accelerate 
the fin either forward or backward, depending on its sign. For 
the above conditions the mechanism is accelerating, i.e., the 
pressure force is much higher than the viscous force. If the rope 
is broken and the ship is allowed to displace, after accelerating 
it tends towards an asymptotic state in which the mean velocity 
is constant and the net thrust (i.e., sum of pressure and viscous 
force) is zero. For the parameters of the figure above, this 
situation is obtained at U „= 5.1 m/s, and the respective forces 
are represented in Fig. 6. 

From the figures above, it can be observed that the 
instantaneous pressure force reaches its maximum value 
twice in each cycle because of the symmetrical undulating 
movement. The maximum thrust takes place at the time that 
the ends of the fin are in the position corresponding to zero 
amplitude. On the contrary, the worst case occurs at the point 
where the edges of the fin are in the position corresponding to 
its maximum amplitude. 

The period-average thrust and efficiency of the three cases 
studied are summarized in Table 1, in which the negative 
sign was used to denote the direction opposite to that of the 
stream. 
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o-—-— Pressure force 
Viscous force 

—o— Total force 
Average total force 


“0 0.2 0.4 0.6 0.8 l 
Time (t/T) [-] 
Fig. 6. Time variation of instantaneous thrust and drag and average total 


force at U „= 5.1 m/s; f = 10 Hz; A, = 0.02 m and à = I m. Configuration 
(a), amplitude increasing from the anterior part to the posterior 


Tab. 1. Period-averaged thrust and efficiency for the three fin undulating 
patterns at U = 5.1 m/s; f = 10 Hz; Aa = 0.02 mand i= 1 m. 


Fin undulating : > 


Configuration (a) 


Configuration (b) 


0 44.4 
1.8 34.5 


(a) Amplitude increasing from the anterior part to the posterior. 

(b) Constant amplitude. 

(c) Amplitude increasing from the anterior part to the mid part 
and decreasing towards the posterior. 


From Table 1, it can be seen that the case with amplitude 
increasing from the anterior part to the mid part and decreasing 
towards the posterior (configuration c) is the worst configuration 
because it corresponds to the lowest efficiency and thrust. The 
case with constant amplitude (configuration b) produces the 
highest thrust and the one with amplitude increasing from the 
anterior part to the posterior (configuration a) produces the 
highest efficiency. This conclusion was expected because the 
travelling waves with increasing amplitude are a typical form 
adopted for swimming motion of live fish in nature. The most 
effective amplitude configuration, case a, will be discussed 
further in the following sections. 


Pressure field 


Fig. 7 shows the pressure distribution in the zone closed to 
the fin. It corresponds to the pattern with amplitude increasing 
from the anterior part to the posterior at U, = 5.1 m/s; f= 10 
Hz; A ax = 0.02 m and à = 1 m. i 

The iso-pressure contours show areas of high and low 
pressures. Behind the fin a high pressure region is created and 
hence a propulsive jet which produces thrust. As discussed 
previously, the peak for pressure and viscous forces do not 
occur simultaneously. The maximum pressure force occurs 
at instants of approximately maximum amplitude at the tail, 
when the pressure gradient attached to the body acts to pull 
the fish upstream, 0.25t/T and 0.75t/T in Fig. 7. On the other 
hand, the minimum viscous force occurs when the undulatory 
wave has a maximum amplitude at approximately the mid- 
point of the fin, 0.5t/T in Fig. 7, when the pressure gradients 
are smaller. 

From Fig. 7 it can also be seen that a staggered array of 
reverse von Karman vortices forms behind the fin, which 
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0.25t/T 


0.5t/T 


0.75UT 


Fig. 7. Pressure field (Pa) at U „= 5.1 m/s; f= 10 Hz; A 
and i = I m. Configuration (a), amplitude increasing 
from the anterior part to the posterior 


= 0.02 im 


max 


decay downstream due to viscous dissipation. These vortices 
are very important because they are closely related to the jet 
stream which generates an opposite force to thrust the body 
forward and their optimal formation would increase efficiency 
noticeable. 


Velocity field 


Fig. 8 shows the velocity field of the pattern with amplitude 
increasing from the anterior part to the posterior at U, = 5.1 m/s; 
f= 10 Hz; A „= 0.02 m and à = 1 m. It corresponds to Fig. 7. 

6.00e+00 2 y Po ne A A x = im 

5.906400 [= = 
5.80e+00 | 

5.70e+00 |- 
5.60e+00 
5.50e+00 
5.40e+00 
5.30e+00 
5.20e+00 
5.10e+00 
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4.90e+00 
4.80e+00 
4.70e+00 
4.60e+00 
4.50e+00 
4.40e+00 
4.30e+00 
4.20e+00 
4.10e+00 
4.00e+00 


Fig. 8. Velocity field (m/s) at U, = 5.1 m/s; f = 10 Hz; A 


=(.02m 


max 


and à = 1 m. Configuration (a), amplitude increasing 
from the anterior part to the posterior 


In this picture, the downstream-directed jet and the von 
Karman vortices are appreciated too. The streamwise velocity 
in this thrust jet is as much as 20% greater than the free-stream 
velocity, so the fish increases the momentum of the fluid in 
the wake. 


Influence of the amplitude, 
frequency and wavelength 


The influence of the maximum amplitude, oscillation 
frequency and wavelength on the thrust and efficiency was 
analyzed and represented in Fig. 9. This figure was carried out 
for the configuration (a), amplitude increasing from the anterior 
part to the posterior and zero free-stream velocity. The other 
two configurations presented similar tendencies. 


Ce-- Efficiency [%] 
O---- Thrust [N] 


Average Thrust [N] 


O---— Efficiency [%] 
O----: Thrust [N] 


Average Thrust [N] 
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Alm] 


O—-- Efficiency [%] 
O----- Thrust [N] 


Average Thrust [N] 


0.25 0.5 0.75 1 1.25 1.5 
à [m] 
Fig. 9. Variation of the period-average thrust and efficiency with amplitude, 


frequency and wavelength at U_= 0 m/s. Configuration (a), amplitude 
increasing from the anterior part to the posterior 


From the figures above it can be seen that, as expected, the 
thrust increases when the oscillation frequency increases. As 
the frequency increases, the efficiency increases rapidly to its 
maximum and then decreases gradually. Other studies realized 
in flapping foils and fishes confirmed that the efficiency is 
maximum in a limited range of frequency. (Triantafyllou, 1995; 
Liu and Kawachi, 1999; Deng et al., 2005). 

In what concerns the maximum amplitude, it can be seen 
that the thrust increases with increasing amplitude. However, 
this it at the expense of even larger power cost required and 
hence less efficiency. On the other hand, at smaller amplitudes 
the efficiency is higher but it is more difficult to generate an 
effective thrust. 

In what concerns the wavelength, as the wavelength 
increases, the average thrust at first increases up to a certain 


maximum value and then decreases, while the efficiency 
increases monotonically. 

From the calculated results, the best option should be 
a combination of thrust and efficiency. It would depend on 
every practical case. For example, military applications require 
high specific power even though the efficiency is not optimal. 
On the contrary, for civil applications the most important issue 
is the efficiency. 


EXPERIMENTAL RESULTS 


A validation test was carried out in order to assess the 
adequacy of the computational model developed here. The 
experimental results were based on tying the ship to a fixed 
point by means ofa rope, Fig. 5b. Once the prototype was fixed, 
the force was measured at different oscillation frequencies 
of the fin. A load cell was used to measure these forces and 
a variable-frequency drive control was used to control the 
oscillation frequency of the electrical motor. 

Numerically, it was simulated with zero free-stream 
velocity. Numerical and experimental results for the average 
thrust force against the frequency are compared in Fig. 10. 
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H 
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Frequency [Hz] 


v Experimentally obtained thrust 
Numerically obtained thrust 


Average Thrust [N] 


6 7 8 9 10 11 12 13 14 15 
Frequency [Hz] 
Fig. 10. Average force against frequency for numerical and experimental 


results. a) amplitude increasing from the anterior part to the posterior; 
b) constant amplitude 


From the figure above, it can be seen that the experimental 
results were in good agreement with numerical ones, which 
means that this CFD model can accurately predict the physical 
phenomena. 


CONCLUSIONS 


In this work, a fish-like marine propulsor was designed and 
the fluid around the mechanism was numerically investigated 
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by solving the three-dimensional incompressible Navier-Stokes 
equations. The hydrodynamics of the fluid flow, pressure 
and velocity fields were obtained. The efficiency and forces 
(which were decomposed into pressure and viscous forces) 
were obtained too. The influence of the amplitude envelope 
was studied, concluding that the configuration with amplitude 
increasing from the anterior part to the posterior is the most 
efficient. Several values of the frequency, wavelength and 
maximum amplitude were analyzed and it was found that the 
efficiency and thrust depend strongly on these parameters. 
Particularly, the thrust increases both with the increase of 
the frequency and the amplitude; the efficiency decreases 
with increasing the amplitude and as the frequency or the 
wavelength increases, the efficiency increases to its maximum 
and then decreases gradually. For this reason, the results 
obtained in this study provide interesting information on the 
performance of fish-like mechanisms and demonstrate that 
a properly design of the fin can be very appropriate for ships 
and underwater vehicles. The optimal propulsive mechanism 
would be a combination of the properly thrust or efficiency, 
depending on every practical application. 

While CFD is a very powerful tool, it can not be applied 
blindly as it has the capability to produce non-physical results 
due to erroneous modeling. For this reason, the numerically 
obtained thrust was compared with experimental results from 
the laboratory mechanism to ensure that the model accurately 
predicts the physical phenomena. 
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NOMENCLATURE 


— Fin length [m] 

— Time [s] 

— Force [N] 

— Amplitude [m] 

— Surface [m?] 

— Unit normal vector [-] 

Velocity [m/s] 

— Period [s] 

— Oscillation frequency [Hz] 

— Pressure [Pa] 

— Gravity [m/s’] 

Empirical constant coefficients [m7/s”] 
— Turbulent kinetic energy [m7/s?] 


FOMBD He eB USHI 
Q 
| | 


= 
m 


— Viscous stress tensor [N/m’] 

— Density [kg/m*] 

— Dynamic viscosity [Pa-s] 

— Wavelength [m] 

Efficiency [-] 

— Kronecker delta [-] 

— Turbulent viscosity [kg/m-s] 

— Viscous stress tensor [N/m?] 

— Dissipation rate of turbulent kinetic energy [m7/s”] 


“At HS PEDA 
| 


Subscripts 

x,y,z — X,Y, Z directions 

i — Cartesian coordinate (i = 1, 2,3 
j — Cartesian coordinate (j = 1, 2, 3 


16 POLISH MARITIME RESEARCH, No 1/2011 


max — Maximum 
p — Pressure 

v — Viscosity 
T — Total 

S — Swimming 
U — Useful 

eo) — Free-stream 
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ABSTRACT 


Added mass is an important and effective dynamic coefficient in accelerating, non uniform motion as a result 
of fluid accelerating around a body. It plays an important role, especially in vessel roll motion, control 
parameters as well as in analyzing the local and global vibration of a vessel and its parts like propellers 
and rudders. In this article, calculating the Added Mass Coefficient has been examined for a sphere, 
ellipsoid, marine propeller and hydrofoil; using numerical Boundary Element Method. Since an Ellipsoid 
and a sphere have simple geometric shapes and the Analytical values of their added mass coefficients 
are available, so that the results of added mass matrix are obtained and evaluated, using the boundary 
element method. Then the added mass matrix is computed in a given geometrical and flow specifications 
for a specific propeller and its results are studied versus experimental results, which it’s current numerical 
data In comparison with other numerical methods has a good conformity with experimental results. The 
most important advantage of the method in determining the added mass matrix coefficients for the surface 
and underwater vessels and the marine propellers is extracting all the added mass coefficients with very 
good Accuracy, while in other numerical methods it is impossible to extract all the coefficients with the 
Desired Accuracy. 


Key words: Added mass coefficient; Boundary element method; Marine propeller; Ellipsoid 


INTRODUCTION 


In fluid mechanics, added mass or virtual mass, is the added 
inertia to the system, owing to that the Increase or decrease in 
the body acceleration should cause the Fluid to move around the 
body in such a way that the object can move through it, because 
the body and the fluid cannot simultaneously occupy this 
physical space. For simplicity this can be modeled in such a way 
that the Volume of fluid moves along the object, While in fact 
all parts of the fluid move with varying degrees. Added mass 
coefficient matrix usually is expressed as a dimensionless value 
which is equal to added mass coefficient divided by the Mass of 
displaced fluid. For example the body density can be considered 
in the body volume. In general, the added mass coefficient is 
a second order tensor that depends on the fluid acceleration 
vector in order to obtain force vector on the body. 

The added mass was considered as the first time in 1776 by 
Dubua who did laboratory studies on spherical pendulum in 
low swings. [2] Precise mathematical equations for the added 
mass of sphere, was obtained by Green (1833) and Stokes 
(1843) respectively. Stokes, also studied on the sphere motions 
ina finite volume of fluid. Later, following the efforts of many 
researchers, the concept of added mass was generalized for 
moving a desired object in different flow regimes. [3] 


THE ADDED MASS COEFFICIENT 
RELATIONS 


The hydrodynamic forces and moments are determined by 
fluid inertia and viscous properties according to the motion of 
a body in real incompressible fluid. Definitely, the other forces 
and moments can be calculated by obtaining one of the inertia 
or viscous terms. Since the fluid can be assumed as an ideal 
one (Inviscid). Therefore the inertia forces and moments can 
be calculated with this estimate. The forces and moments of 
inertia can be introduced as the body added mass terms, using 
this method. 

Assume that the body with surface S is moving in an infinite 
homogeneous ideal Fluid. Consider two systems of coordinates: 
the stationary one (We denote it by XYZ) and the coordinate 
system moving together with the body (We denote it by Oxyz). 
The vortex-free condition implies the existence of a potential 
o(X, Y, Z, t), assuming that the moment of these two systems 
of coordinate are specified. So the components of the fluid 
velocity are defined as follows: 


VS" > Yonah i oe (1) 
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Using the continuity and incompressibility of fluid and due 
to the vortex-free condition as mentioned above the Laplace 
equation is implied as follows: 


p 2e Fo 
Xx? oY oZ 
The boundary conditions for Eq. (2) look as follows: 


=0 ; Vþ=0 2) 


1. The watertight condition, valid on the surface S: 
op 
on|s 


| =U,) is a projection of fluid velocity on the 


=u (3) 


n 


Where ce 


(external) direction n, which is normal to surface S; u, is the 
velocity projection of a point of the body, normal to surface n. 


2. Stationary condition at infinity: 


ime” = eed = im =() (4) 
r>» OX r>% OY roo QZ 


Where r? = X? + Y? + Z? (r is the distance from the origin to 
a fluid point). The function vanishes at infinity as 1/r’, whereas 
its first-order coordinate derivatives vanish as 1/r? [4, 5] 

From the formulation of the problem (1-4) it is seen that 
the function @ depends on time t via the boundary condition 
of right-hand side equation (3), so the boundary condition is 
investigated specifically. The origin O is chosen to coincide 
with an arbitrary point of the body 

The velocity of the point O is Denoted by Up (components 
Up Upy Up, are projections of the vector u, on coordinate axes 
attached to the body). 

The angular velocity of the body with respect to the point O 
is also donated by @ (components of @ in the same coordinate 
system are denoted by o,, @,, @,). Thus, the velocity of an 
arbitrary point of the body, including any point of its surface 
S, is determined by the following equation: 


u=u,+@xT (5) 
Where f is the vector, determining the position of the point. 
The components of Eq. (5) look as follows: 
U, =U, +@,Z—- O,Y 
u, =U, +0,X — O,Z (6) 
u, =Uy,+@,Y — O,X 
On the surface S we have: 
u,=u,cos(n, X)+u,cos(n, y)+u,cos(n, Z) (7) 
Writing: a = cos(n, x); B = cos(n, y); y = cos(n, z) and 


substituting (6) into (7), the following form of the boundary 
condition is obtained (3): 


al = u= Up, ot UB oe Uo,Y T oO, (yy T zB)+ 
On|. (8) 


+0, (za —xy)+ @, (xB - yo) 


In the formula (8) the variables a, B, y, yy — zB, za — xy, 
xB — ya are determined only based on the shape of the body 
surface. The body motion and the dynamics of the flow 
are determined by the functions u, Up Uy, ®,, ®, ©, The 
potential ọ is represented as follows, due to the Linearity of 
the problem. 


P=Uy, 9, T Up, Py T Uoz93 + O94 + OyOs+ QRO (9) 
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In the above equation f, i= 1, 2, 3 are the flow potentials 
corresponding to the body moving along the axes x, y, z at unit 
linear velocities, respectively. On the other hand 6, i = 4, 5, 6 
are potentials, corresponding to rotation of the body around the 
same axes at unit angular velocities respectively. It is seen that 
the body motion in an ideal infinite fluid causes to solve six 
problems. The first problem can be formulated as follows: the 
solution of the Laplace equation Ag, = 0 will be found, using 
the following boundary conditions: 


0, 


= 10 
aoe (10) 
Mi Ns GR, as r— œ (11) 
ox OY Z 
a 


The function —1 r — =q is given on the surface S and The last 


(sixth) problem, taking (8) into account. To solve the problem, 
considering the Laplace equation Ad, = 0 using the following 


boundary conditions: 2h =xß-ya which is given at the 


surface S, so: 
Ae | Abe | Be 


oX OY OZ 


It can be seen that the function Q; i= 1, 2,... 6 does not 
depend on u, and @,. These functions are determined only by 
the shape of the surface S of the body and Choice of coordinate 
system Oxyz attached to the body. 


KINETIC ENERGY OF THE FLUID 


>0 as r> o 


In order to to obtain fluid kinetic energy, considering the 
surface S and the sphere È of radius a including the body with 
surrounding Fluid, the following integral is defined. 


retofffyav- 
(12) 


=5p Í f Í, (2) (2) (2) dxdydz 


Where p, is the fluid Density 
Using Green’s transformation for two functions (@,, 0,): 


if (> h , A, 2 , A, Ob, e lix ids 
Ox ax Oy Oy Oz Oz 
(13) 


=-Í La das- f ffo, a oo ue 2 © asaya 


And considering , = 6, = 6, Ad = 0 and putting words in 
the above equation the following expression is gotten: 


Tean [o2a dS- 5 She Sas (14) 
2478 ðn 
The Second term in the right hand side = Eq. (14) tends 
to zero as a — œ according to the asymptotic behavior of o 
and its first derivatives while r — œ. Therefore, the following 
formula is obtained for the total kinetic energy of the fluid 
around the surface S: 


a 
e ffo zas (15) 


Substituting the expression (9) in to (15) and writing: 
Up, =Uj, Ugy =Uy, Ugy =U, 


O,= Uy, O,= Us, @, = Us 


The expression (16) : ae! obtained: 


T= SD u, (16) 
Where i, is = 
ad, 
hiac=—P [jnas (17) 


à, is called the added mass of the body. 


According to Eq. (17) it is found that the added mass, do 
not depend on the body velocity and Kinematics of the motion 
but is only a function of the body surface geometry. 

Applying Green’s formula to functions , and 6, in the 
volume V between the surfaces ÈX and S then: 


Mh, (0, Ad, —0,A,)dV = 
(18) 


7 {i Me n as- ff (4, Bo, 2 hs 


Since Ne =Ao, = 0, so that the left-hand side of the equation 
equals zero, and the first term of the right-hand side of the 
equation tends to zero at a — œ. Therefore, the condition for 
the infinite fluid surrounding the body is as follows: 


E A Bas aw 


Since i, = à; Therefore, only 21 values out of 36 A, (i, k= 
= 1,2... 6) are independent. 

This part, in fact proved that the added mass matrix is 
symmetric. 

To obtain the added mass coefficient matrix using force 
view, ultimately the same formula is achieved. [6] 


CALCULATING THE ADDED MASS 
COEFFICIENT USING THE BOUNDARY 
ELEMENT METHOD 


Boundary element method is based on the Green’s theory. 
Based on the theory solving the fluid flow around the body can 
be transferred from the volume solving domain around the body 
to the body boundary that is performed on the produced surface 
elements in three-dimensional modes. The boundary element 
method is used to obtain all six unit potentials corresponding to 
the six free degrees. The added mass tensor can be computed, 
obtaining the unit potential. The forces and momentum also 
can be achieved, having the added mass. [7] 


fV.udo= funds 


Q F 


THE BOUNDARY ELEMENT METHOD 
TO COMPUTE THE ADDED MASS 
COEFFICIENT 


(20) 


Consider the Closed area Q, the boundary S and unit vector 
n normal to S (Fig. 1). The boundary S is composed of the 
body wetted surface SB, the wake surface SW and the external 
control surface S, including surface area SB, SW. 

The Q encountered the Uniform inflow velocity Ų, 
assuming the flow is Incompressible, Inviscid and Irrotational 
in the domain Q, in the upstream flow. 


S << 


Fig. 1. Application of Green theory for hydrofoils 


With these assumptions, the flow field around the body can 
be identified using Perturbation Velocity potential þ, which 
satisfies the Laplace equation: 


Vo =0 (21) 


To solving the problem, the boundary conditions are 
considered as follows: 


Kinematic Boundary Condition 


Based on the Kinematic boundary condition, the Flow 
velocity normal to the body surface must be zero. 


—=0 > cal =- V.n 
on on 
Where ñ, is the unit vector normal to outside of the boundary 
[8] and is changed as follows, considering the unit potential 
Kinematic of Boundary Condition: 


(22) 


ee 

0; _ i]j=133 (23) 
j=436 

on Pie (xn) ? 


i=l>3 


Where, r is a vertical distance from a fixed origin point 
[12]. 


Kinematic And Dynamic Conditions In Wake 


Based on this condition, there is no velocity jump of flow 
on the wake vortex, however there is velocity jump of potential 
on the surface which is equal to circulation I around the fin. 
The mathematical relations of flow velocity and potential are 
expressed as follows: 

=o°-¢ =T 


(Apa s, 
a- yp Ay" = 


Where B and F a ai and = sides of the propeller 
[8]. 


KUTTA CONDITION IN THE TRAILING 
EDGE 


Hess and Smith have shown that the flow passing along 
a thin and non lifting body can be described by the distribution 
singularities of source. But to describe the flow passing along 
lifting bodies, modeled by the circulation distribution on the 
surface, a boundary condition for the Trailing Edge of the body 
must be satisfied, expressing that the velocity at the trailing 
edge should be limited and unique. [9, 10] 


[Vol <0 (25) 
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Surface At Infinity 


At the boundary surface at infinity S_, the perturbation 
velocity due to the body surface must tend to zero. 


Vo —> 0 as S, —> 0 (26) 


Using Green’s theory, general answers of the Laplace 
equation for each point P can be expressed as the following 
integral equation: 


An o(P)= jpa oe )ds+ 


(27) 
o>. 1 1 
“Shen, Rea ce s+ [fabs aop" 


In the above equation ô/ðn is determined, using the 
boundary condition of Eq (23). 1/R(p, q), is the distance 
between p and q, so that the only remaining unknown for 
solving Eq (27) is . [8] 


EQUATIONS DISCRETIZATION 


In order to solve Eq (27), the equation form must be discrete 
and the body boundary surface shall be made of tetragonal 
elements. 

Fractional equations lead to linear algebraic equation system 
for p and i as follows: 

Ne Nw 


20; = >Pse J+ D Wi (Ad), + 


Esch 


Where D;, Sj, and W; are Dipole coefficients and source 
for j element that acts on the arithmetic point i. The coefficients 
are defined as follows: 


(28) 
i=12:.,N 


tot 


DeU T )ds,], on S, 
w => Z S Jds; on Sy ©) 
S; -SLf ash on Sp 


To obtain these coefficients the numerical solution of the 
integrals is used. Calculating these coefficients, the following 
matrix equation is formed 

Solving the equations using the Gauss - Sydl method is done 
as a repetitive solving, so that Q; are obtained. [11] 


[D]lo] = [S] + [W][A9] (30) 


Where: 

[D] = dipole Potential effects of induction coefficient matrix 
[S] = source Potential effects of induction coefficient matrix 
[W]= dipole Potential effects of induction coefficient matrix 


For lifting and non-lifting bodies, the above computational 
method can be used, but the difference is in the use of Kutta 
boundary condition which there is no need to consider it in 
non-lifting bodies, so that the results are obtained directly 
regarding the kinematic boundary condition. 
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Although the Kutta condition is considered to analyze 
the propeller and the hydrofoil, the obtaining results have no 
significant difference with the ones which are obtained without 
considering the Kutta condition. So, regarding this issue it can 
be expressed that only the D and S out of the D, S and W are 
important for calculating the added mass matrix. 


k ð -1 
D; POR on Ss 
k=l § j ij 
7 (31) 
k 
1 
S; = LUIGI on Di 


Which in Reference 12, the integrals are expressed as 
follows that in fact are the same general term in Eq (27) and with 
extracting the D and S, the formulas in Eq (31) is obtained. 


1 Lg Sl 
IF #hs-5 i ôn 


K#i j=!>6 


3S, 


ofl “ofi a 
| 02 (2) = 21) dS, + 27;_« 
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Using The Boundary Element Method 
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matrix 


If the above integrals are formulated based on matrix 
equations, the following relation is obtained: 


[A diod] = [B] 
So that i = K and A x = 21. 
The added mass matrix values with solving matrix equations 
for all six potential unit functions can be defined as follows: 


Hix =P fonds, 
Sg 


If the above equation is written in the matrix form, which 
is the general form of added mass matrix. [12]: 


Added Mass Matrix = 


(33) 


(34) 


N 4 . N . . 
Doim,58, .. D1 g)(rx 038, 
i=l i=l 


N N , 
S pn 8S, ... A dr Rn) iS, 
i=l i=l 


p“ i ae 
So 6n58S; . $0 xn) 8S; 
i=l i=l 
RESULTS 


Computing The Added Mass Coefficient 
Of Sphere 


The added mass coefficient for a sphere of radius 0.25 m 
has been investigated in this part. However, according to the 


analytical relations in the mass matrix of sphere, only a,,,a,,and 
a,, have values and the others are zero. The three coefficients 
are equal with each other and should be equal to half of the 
sphere mass. In order to evaluate the Mesh independency, the 
added mass matrix were investigated for three cases of surface 
element (2400, 3360 and 4800 elements), using the boundary 
element numerical code and the obtained results and errors 
were checked. 

The obtained values for the added mass of the sphere using 
numerical analysis of the boundary element are as follows: 


Added Mass Matrix = 
0.032752 0 0 0 0 0] 
0 0.032753 0 0 0 0 
0 0 0.032753 0 0 0 
“| 0 0 0 000 
0 0 0 0 0 0 
0 0 0 000 
For a sphere with radius of 0.25 m, the analytical values of 
aip 4,5, a, are equal to 0.0327245 


As it was mentioned, only a,,, a,,, a,, have non-zero 
value, which is equal to half of the sphere mass and the 
other elements of the added mass matrix are zero. So it can 
be seen that the obtained values of the numerical boundary 
element method have a high conformity with the analytical 
data and the error rate, according to the table (3) is the 10E-2 
order. [6] 


Tab. 1. Geometric and flow data to calculate the added mass coefficient of sphere 


Radius = 0.25 
Volume (Numerical) = 0.065449 
Volume (Exact) = 0.0654498 


Fluid = Water 
Density = 1000 kg/m? 


Tab. 2. Showing the produced geometry with surface elements to check the mesh independency condition 


4800 surface elements 


3360 surface elements 


3360 surface elements 


Tab. 3. The compared Results of the added mass numerical and analytical data of the sphere and the Error percentage Due to the elements numbers 


Number of | Numerical result for 
Element a,, (code result) 


0.032814 


Analytical 
value of a,, 


0.032725 


Difference between 
Numerical and Analytical result 


8.95E-05 


% Error 
411(Num)~ 4 11(Anal) 


)x100 
a 11(Num) 


0.272749 


0.032771 


0.032725 


4.65E-05 


0.141894 


0.032752 


0.032725 


2.75E-05 


0.083964 
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Tab. 4. Geometric and flow data of computing the added mass coefficient ellipsoid 


Volume(Numerical) = 0.167376 
Fluid = Water 
Density = 1000 kg/m* 


4600 surface elements considering the mesh independency 


Computing The Added Mass Of Ellipsoid 


Continue the contex, the coefficient of added mass for an 
ellipsoid has been investigated. 

Since the added mass coefficients of ellipsoid as well as 
the sphere are available in references, so that to validate the 
computing results of the boundary element code, this analysis 
is given. The numerical analysis with the mesh independency 
condition was investigated and the results of the surface element 
with the number of 4600 had very little error in comparison 
with the existing data. 

The added mass matrix obtained by numerical analysis for 
an ellipsoid with the above specification is as follows: 


Added Mass Matrix = 
0.009025 0 0 0 0 0 
0 0.143401 0 0 0 0 
| o0 0 0.1434010 0 0 
| 0 0 0 0 0 0 
0 0 0 00.030467 0 
0 0 0 0 0 0.030467 


Considering that the numerical analysis results of the added 
mass were dimensional and the results in the references were 
dimensionless, then to validate the results, the data obtained 
from the numerical analysis were transformed to dimensionless 
one and compared with the available data in those references. 


To make the added mass matrix of ellipsoid dimensionless a, ,, 


a,, and a, shall be divided on (42/3)ab? and also a Were and a 
must be divided on (42/15)ab? (a? + b°). 
Due to the symmetry of the ellipsoid, then a,,=a,,, a,,=a,, 


and a,, = 0. 
According to data obtained from the Diagram and making 
them dimensional, then: 


4 
a =K; xab’ 


4n > 
Ase = ay = Kay an 


aig Hane Re xaa’ +b’) 


Itis found from comparing the boundary element numerical 
analysis results with the data available in references for added 
mass coefficients of ellipsoid that there is a good conformity 
between the numerical data and the ones in the references 
[6, 12]. 


Tab. 5. The added mass coefficient Diagram of an ellipsoid [6] 


Non Dimentional value 
kıı = 0.05384 
kz = 0.85384 
ks5 = 0.89538 


4/3nab* = 0.167551 
4/15mab*(a°+b”) = 0.034850 


Element numbers were 4800, considering the mesh independency 
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Tab. 6. the compared Results of the added mass numerical and analytical data of the ellipsoid and the Error percentage Due to the elements numbers 


% Error 
411(Num)~ 2 11(Anal) 


Difference between 
Numerical and Analytical result 


Numerical result 
(code result) 


Analytical 
result 


)x 100 
a 11(Num) 


0.044986 
0.236581 
-2.41898 


4.06E-06 
0.000339 
-0.00074 


0.009025 0.00902094 
0.143401 0.14306174 


A 0.030467 0.03120399 


Computing The Added Mass Matrix For A Sample Propeller 


According to the data and results obtained for the added mass matrix of the sphere and ellipsoid it can be found that there 
is a very good accuracy between the data obtained from the analytical and the numerical solution, suggestive of the validity of 
the boundary element method code in calculating the added mass matrix.The following text is about computing the added mass 
matrix for a sample propeller which its experimental and geometrical data is on hand and the results have been validated with 
the experimental data.The added mass matrix, considering what was mentioned in the text do not depends on the body velocity 
and is only a function of the body geometry. 

First of all, the existing Geometric data of the propeller was considered based on the standards and the propeller geometry 
modeled with a very high accuracy. Then the propeller surface and the hub were analyzed according to the flow conditions. The 
propeller analysis was performed in several different modes of surface element numbers so that the most appropriate numbers 
of the elements is achieved in a way that the obtained answers do not depend on the element numbers [14]. 


Tab. 7. Geometric and flow data related to the propeller and the produced geometry for boundary element numerical analysis [13] 


Type FPP 
No. of Blades, Z 5 


Expanded Area Ratio 0.595 
Advance Ratio, J 0.4489 
Diameter, D (m) 0.25 
Material ALBC3 


Young's Modulus, E (GPa) =I 


Shear Modulus, G (GPa) 
Material Density, pp (kg/m3) w 
Water Density, | Water Density, pw (kg/m3) | | Water Density, pw (kg/m3) | 


Revolution Speed, n (rpm) — 4 


Dimensional Added Mass Matrix Obtained For The Propeller 


2.176192 -0.240440 -0.686280 -0.670630 0.015638 -0.196440 
-0.240440 0.042111 0.080637 0.006256 0.000156 0.018995 
Added Wiest -0.686280 0.080637 0.251128 0.019366 -0.005240 0.059255 
-0.058640 0.006256 0.019366 0.001617 -0.000450 0.005282 
0.015638 0.000156 -0.005240 -0.000450 0.000505 -0.001630 
-0.196440 0.018995 0.059255 0.005282 -0.001630 0.018315 
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Tab. 8. Compared results of the added mass numerical and analytical data of the propeller and the Error percentage Due to the element numbers 


Numerical result Diffrence between % Error 


Number of x Exprimental : ; a11 (Num) 7 2 11(Anal) 
Fiement Matrix Array A, Result Numerical and Exprimental um na )x100 


(Addedprop code) result a 11(Num) 


0.6242060 0.6860000 -0.0617940 -9.0078717 
0.6382161 0.6860000 -0.0477839 -6.9655831 


0.6451033 0.6860000 -0.0408967 -5.9616181 


0.6510313 0.6860000 -0.0349687 -5.0974781 


0.6631462 0.6860000 -0.0228538 -3.3314577 


0.6691640 0.6860000 -0.0168360 -2.4542274 


0.6706359 0.6860000 -0.0153641 -2.2396647 


Investigating The Mesh Independency 


To investigate that the obtained data for the propeller is independent of the produced surface elements of the propeller and the 
hub surface, the added mass matrix was obtained in seven different numbers of surface elements and the results were checked, 
as it is identified from the diagram (figure 2) the results do not have any significant change after 3960 elements. 


rrr rri 


ra R | 


Added Mass Coefficient A13 
Ò 
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1 1 L L 
2000 3000 4000 5000 
Number of Element 


Fig. 2. Changes in the mass added coefficient based on the element numbers at the propeller surface 


Comparing The Added Mass Matrix Of The Propeller With The Obtained Data From Other 
Numerical And Experimental Methods 


Tab. 9. Comparing the data, obtained from the numerical solution, experimental and formulated data 


Added Mass 
(kg) 


Experiment 0.6860000 


Estimation Method 


AddedProp (Present method) 0.6706359 


Prodas 0.6690000 
Schwanecke 0.8240000 


Shuster 0.9660000 
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Computing The Added Mass Matrix Of The Hydrofoil (NACA0012 Section) 


Tab. 10. Geometrical data and the element numbers of calculating the added mass coefficient of hydrofoil 


Length of Span = 1.0 
Length of Min Chord = 0.9 
Length of Max Chord = 0.9 
Number of Element = 2632 
Section profile: NACA - 0012 


Gradual Increase In The Span Length To Validate The Results 


To calculate the added mass matrix of the hydrofoil with NACA0012 section, first of all the mesh independency was 
investigated in which the obtained data from 2632 elements had very good accuracy, then to validate the obtained data, the span 
length was increased gradually from 0.2 to 1 and thus the added mass matrix was computed. Afterward the added mass matrix 
element changes were plotted based on the span length increase. All of the coefficients are increased with linear and non-linear 
curves of the 2 and 3 order, which is completely ascending. In fact, by increasing the Span length, all coefficients were reasonably 
increased as expected. Of course, investigating each of the diagrams whether they should be linear or nonlinear, of the 2 orders 
or 3 orders are denied due to that they are Proportional to Dynamic and Simulation topics. 


—G— All 
0.15 | — A— Al2 
| —--7-—- A15 


= ALE 

f = 4 — 422 
=—--Q-—- A33 

r —-O— ade 


L — = — A55 
+ A66 


Added Mass Coefficient 


06 
Lenght of Span 
Fig. 3. The added mass coefficient changes of hydrofoil based on the gradually increase along the Span length 


The Added Mass Matrix Of Hydrofoil 


The added mass matrix of hydrofoil with NACA0012 section for the Span length of 1 m and Chord length of 9/0 m is as 
follows: 


-13.5806 89.37951 58.75352 
-63.4481 11.35911 44.82613 
0 0 0 
40.92427 -7.51047 -23.5442 
-7.51047 56.3155 21.46682 
-23.5442 21.46682 22.75078 


159.3152 31.63824 
31.63824 122.4136 
Added Mass Matrix = 0 0 

-13.5806 -63.4481 
89.37951 11.35911 
58.75352 44.82613 


O O O O O OO 
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CONCLUSION 


e The added mass coefficient is one of the most important 
dynamic coefficients in solving dynamic equations and 
simulating the vessels. Calculating the Add mass coefficient 
especially the bodies whose shapes are not simple is 
difficult. Added mass coefficient is extractable using 
empirical formulas, the simple analytical relations and 
numerical methods such as Strip theory and lab tests. 

e But each of these methods can be used in a limited range. 
For example, the analysis relations can be used only for 
simple geometries and lab tests to extract the added mass 
are costly and obtaining all matrix coefficients is very 
difficult. Moreover, the numerical methods cannot model 
the full three-dimensional body and have a lot of errors 
that ultimately causing a large deviation in calculating 
the coefficients. Meanwhile since the boundary element 
method solves the problem on surface boundary elements 
of the body, it needs less elements and time to solve and 
can obtain the added mass coefficients owing to modeling 
the full three-dimensional body and also extracting the 
potential values with high accuracy. 

e As mentioned in the article to investigate the mesh 
independency, the added mass matrix of an ellipsoid and 
a sphere of specific geometry were extracted, using the 
boundary element method and the results were validated 
with the analytical data. 

e Afterward the added mass matrix for a certain propeller was 
computed which its results had a very good conformity with 
the experimental ones. However due to lack of access to 
distribution Cord of the propeller, the standard distribution 
of Series B for a propeller with specific section (NACAa66) 
was used, considering that the volume of the produced 
propeller geometry was checked with the volume of the 
examined propeller in a way that the geometry would have 
a high conformity with the existing data so that if the exact 
geometry is on hand, the error rate with be lower than 10E-2 
order. 

e To put it inanut shell, extracting a complete and dimensional 
add mass matrix with any desired geometry and the accuracy 
of obtained data is the difference between this method and 
boundary element code with other numerical methods. 
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A method for evaluating durability 
of porous slide bearings 


Piotr Bzura, Ph. D. 
Gdansk University of Technology 


ABSTRACT 


The article presents the structures of porous slide bearing bushes and greasy lubricants, 
and the similarity criterion resulting from these structures. A physical model is proposed 
which describes the operation of the porous slide bearing bush taking into account, via 
analogy, the phenomena taking place in the greasy lubricant nodes instead of those observed 
in the bush structures. 
Included are the results of experimental tests of three greasy lubricants which allow 
transferring the results recorded in model experiments to real objects. 


Key words: structure; durability; strength of bonds 


INTRODUCTION 


The published data on operating parameters of Diesel 
engines used in shipbuilding reveal a remarkably slowing 
increase, or even decrease, of the consumption of lubricating 
media [17, 18]. This effect results, on the one hand, from ship 
owners’ permanent search for the ways to reduce the cost 
of operation of the internal combustion engines, and, on the 
other hand, from increasing pressure toward the reduction 
of environment pollution caused by those engines [5]. The 
rationalisation of the use of lubricating media requires 
delivering them to tribological nodes in optimal volumes and 
precisely defined time intervals [11, 12]. 

This need for rationalisation of lubrication is the reason 
why slide bearings with “limited supervision” are expected 
to be used in future in marine piston engines. To this group 
belong, for instance, the bearings with porous bushes. 
A distinguishing feature of these bearings is their durability, 
which can be evaluated in laboratory tests using the here 
presented method. 


DURABILITY OF POROUS SLIDE 
BEARINGS 


The porous slide bearings belong to the group of bearings 
which do not need lubricating service, and therefore they are 
frequently called the slide bearings with “limited supervision”. 
The lubricant is not delivered to them from outside, but it is 


collected in bush wormholes and is to be periodically refilled. 
Taking into account the lubricating supervision, a distinction 
should be made between the porous slide bearings and those 
with conventional lubrication - Fig. 1 [11, 12]. 


Bearing 
slide 


Fig. 1. Scheme illustrating the classification of slide bearings due to 
lubricating supervision: a) slide bearings with conventional lubrication, 
b) porous slide bearings (the lubricating system and the lubricant are within 
the journal-bush area) 


The bushes in the porous slide bearings (Fig. 2) resemble 
a sponge, due to the presence of wormholes having the form 
of capillary ducts [4]. After saturating the porous bush with 
the lubricant using the vacuum method and connecting it with 
the journal, the self-lubricating porous slide bearing is ready 
for operation [9, 11, 12] 
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3 Detail "A" 


Fig. 2. Scheme of porous slide bearings : 1) journal, 2) bush, 3) lubricant 
placed in the bush material, 4) dead end, 5) capillary ducts open at two 
ends, 6) capillary duct open at one end 


The self-lubricity of the porous bearings [3,10,11,12,14,15] 
is realised using the oil flow from the porous wall, which is 
mainly forced by the difference in thermal expansion of the 
lubricating oil and the metal bush structure (Fig. 3). Initially, 
when the bearing is cold, some role is also played by elastic 
strains of the bush material which, pressed by the journal, 
releases the oil from the wormholes (per analogy to a sponge). 
Then the journal moving at the presence of the oil generates the 
hydrodynamic pressure in the oil clearance. In the high pressure 
area the oil is pressed back into the porous wall, then it flows 
through the wormholes and flows out to the oil clearance at 
the low pressure area. 

During the operation, the porous slide bearings lose their 
durability, following the decrease of their lubricating ability. All 
this is caused by oil evaporation and leakage off the bearing, 
due to excessively high temperatures of bearing operation 
and the ageing of the oil. Therefore the area of application of 
the porous slide bearing is limited by its load capacity, i.e. the 
journal load, at which the porous slide bearing can still work 
safely which is affected by [12]: 

— the reduced amount of oil and rate of its circulation, 
— the decreased hydrodynamic uplift pressure. 


Experimental tests are in progress now which are oriented 
on evaluating permissible journal loads in porous slide bearings 
for given sliding speeds at which the slide bearing can still work 
safely [9, 11, 12, 14]. Unfortunately, in these tests different 
criteria of load limits for given speeds are adopted, which 
reduces the space for data comparison and reliability. 

In this context, to provide opportunities for comparing 
journal load limits in different porous slide bearings, a simple 
physical model is proposed in the form of a four-sphere 
apparatus T-02. The frictional association in this machine 
is created using four spheres, each of 12.7 mm in diameter, 
made of LH15 steel in the accuracy class 16 according to the 
PN-83/M-86452 standard. Three spheres are placed in the 
lower holder-container having the form of a small cup, which 
is poured with the greasy lubricant. The fourth sphere, which 


I 


models the journal in the porous slide bearing, is fixed in the 
upper holder. The spheres in the cup are pressed towards the 
sphere fixed in the upper holder using a special lever. The 
parameters controlled during the tests include the journal load 
and rotational speed, and the greasy lubricant temperature 
[19] — Fig. 4. 


Fig. 4. Kinetic scheme of the four-sphere apparatus: I) cover which fixes 
the lower spheres, 2) holder of upper sphere, 3) upper sphere (rotating), 
4) lower spheres (stationary), 5) container with the tested lubricant, 

6) prism, 7) lever, 8) weights, 9) tested greasy lubricant [19] 


DURABILITY OF GREASY LUBRICANTS 


The greasy lubricants [16] are semi-liquid plastic lubricants 
having the colloidal structure which consist of the base oil (the 
mineral oil in general, but in exceptional cases also the synthetic 
oil) and the thickening agent. They are used wherever frequent 
lubricating oil refilling is extremely difficult or impossible. 
They execute all functions of liquid lubricants — except cooling 
and cleaning of friction surfaces. Like for the porous bearing 
bushes, which are the systems consisting of a metal skeleton 
filled with the lubricating substance, the structure of the greasy 
lubricants (composed of the base oil and the thickener, Fig. 5) 
resembles the sponge, and this fact was assumed as the required 
similarity criterion for creating the physical model. 

In the greasy lubricant the percent share of the thickener 
is low, unlike the porous bushes in which the bush skeleton 
plays the role of “the thickener”. Frequently, to strengthen the 
structure of the lubricant and make it more resistant to the action 
of external agents, extra additives are used [16] which modify 
the cohesion forces between the base oil and the thickener in 
the required range. The bigger the cohesion forces, the more 
coherent is the body and the more difficult is to stratify it. 
Therefore the strength of the bonds can be considered the 
property determining durability of the greasy lubricants. 


Fig. 3. Scheme of porous bearing lubrication: I) start-up, II) thermal outflow, II) equilibrium, IV) standstill [12] 
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Fig. 5. Scheme of the structure of greasy lubricants: 
1) thickener, 2) base oil, 3) cohesion forces 


The parameter which determines the durability of 
a conventional slide bearing (Fig. 1a) is the lubricating ability 
of the oil, i.e. the property of the lubricating substance which 
characterises its ability to preserve lubricating features in 
mitigated solid friction conditions. This way it determines oil 
ability to create a stable boundary layer being the consequence 
of the adsorption of fluid particles on the solid body surfaces. 
The measure of the lubricating ability of the oils is the durability 
of the boundary layer, i.e. the durability of bonds between the 
lubricating substance (adsorbent) and the base [1, 2, 6, 8, 16]. 

In this context it seems logical to relate the durability of the 
bonds with the lubricating ability of the greasy lubricants, which 
defines the minimal pressure required for stopping the action 
of the cohesion forces between the thickener and the base oil. 
Therefore the lubricating ability of the greasy lubricants can be 
considered a parameter determining the durability of the slide 
bearing with “limited supervision” (Fig. 1b). 


PHYSICAL MODEL OF A POROUS SLIDE 
BEARING 


The structure of the porous bearing bush is characterised 
by the set of technical and operating characteristics assumed 
when designing and manufacturing the bearing. These 
characteristics determine the level of bearing applicability 
to realise assumed tasks in given operating conditions. Each 
structure has a separate set of structure parameters which map 
the technical condition. 


U= {u} fori=1,,..m (1) 


Change in the technical condition of the porous slide 
bearing is, among other factors, a function of the resistance of 
the porous bush saturated with the lubricant to ageing. During 
bearing operation, the effects of excitation forces acting in the 
micro- and macro- environment of the bush cumulate with 
time, which leads to the leakage flow of the lubricant, and 
bearing seizure as a possible further consequence. Therefore 
an important parameter which determines the durability of 
the porous slide bearing is structural strength of the bushes 
saturated with the lubricant. 

The durability of the greasy lubricant is, in turn, the feature 
which determines the strength H of the bonds between the 
thickener and the base oil (Fig. 5). Consequently, the strength 
of bonds of the greasy lubricant can be indentified with the 
structural strength of the porous slide bearing bush. 

The strength of bonds of the greasy lubricant can be evaluated 
by analysing the diagram of the friction torque (Fig. 6) as the 
time function, created using the results of the investigations of 
the physical model which was the friction pair of the apparatus 
T-02 (Fig. 4) immersed in the greasy lubricant. 
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Fig. 6. Diagram of greasy lubricant friction torque:C,) curve which 
characterises the greasy lubricant revealing small strength of bonds H „ 
C,) curve which characterises the greasy lubricant revealing large strength 
of bonds H, A) time when the boundary layer is broken and reconstructed, 
B) time when the boundary layer is broken and not reconstructed, 

1) upper sphere of the friction pair, 2) lower sphere of the friction pair, 
3) greasy lubricant base oil, 4) friction pair seizures 


The shapes of curves C, and C, depend on the type of the 
greasy lubricant. It results from their nature that the limit of 
strength of bonds between the thickener and the base oil is 
preserved as long as the lubricant resistance to breaking is 
larger or equal to the destroying actions (stage I). Once the 
bonds are broken, the leakage of the base oil takes place. If 
the base oil reveals high lubricating ability, characteristic for 
lubricants [1, 2] and, simultaneously, it has low strength of 
bonds of the greasy lubricant H, (symbol A in Fig. 6), all this 
produces stable bonds between the base oil and the foundation. 
In the opposite case, when the base oil reveals low lubricating 
ability concerning the lubricating oil and, at the same time, 
high strength of bonds of the greasy lubricant H, (symbol B in 
Fig. 6) the boundary layer will not appear and the seizure of 
the bearing will take place. 

That is why the porous slide bearings for which the physical 
model is the friction pair in which the seizure is observed 
(curve C,), are characterised by larger maximal load which 
they can carry. 


EXPERIMENTAL TESTS 


The measurements were performed on the friction pair of 
the T-02 apparatus with spheres immersed in a greasy lubricant 
(Fig. 4). The following lubricants were tested: 

— the graphite grease created on the base of the mineral oil 
with 10% of natural graphite, the minimum, thickened with 
the lime and used for lubricating heavily loaded friction 
pairs in temperatures ranging between -20 and +50°C; 

— the LT43 grease created on the basis of the mineral oil 
thickened with lithium. It contains refining additives, 
especially those acting towards anticorrosion, anti-oxidation 
and improvement of lubricating abilities. It is used for 
lubricating heavily loaded friction pairs in temperatures 
ranging between -20 + +130°C; 

— the grease being a composition of the graphite grease and 
the LT43 grease. 


Due to the flow limits of the greasy lubricants, all tests were 
performed at the temperature of 100°C. Two test procedures 
were executed: 

e determining the strength of bonds H [Nm] at continuous 
load increase for the rotational speed of the spindle equal 
to 500 rev/min, 
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e determining the strength of the structure H [Nm] at 
continuous load increase for the rotational speed of the 
_ spindle equal to 1440 rev/min. 


Friction torque [Nm] 


Applied load [N] 


6 


_8 10 12 14 16 18 
lime tests [s] 

Fig. 7. Joint diagram of greasy lubricants 
for spindle rotational speed n = 500 rev/min 


Applied load [N] 


Friction torque [Nm] 


8 10 12 14 16 18 
Time tests [s] 


Fig. 8. Joint diagram of greasy lubricants 
for spindle rotational speed n = 1440 rev/min. 


The obtained values show that: 

e the smaller the rotational speed of the journal, the larger the 
load P[N] needed for breaking the boundary layer created 
on the spheres, 

e the graphite grease with the thickener having calcium 
compounds reveals very small strength of bonds, and the 
boundary layer between the spheres and the greasy lubricant 
in the friction pair was broken at low seizing load; 

e the LT43 grease with the thickener having lithium 
compounds reveals high strength of bonds and breaking the 
boundary layer took place at higher value of seizing load. 

e the composition of greases consisting in 50% of the graphite 
grease and 50% of the LT43 grease reveals the strength of 
bonds comparable with the graphite grease. 


When the boundary layer is broken, the thickener and the 
base oil separate from each other (points A and B in Fig.6) 
which provides opportunities for evaluating the strength H of 


bonds of the greasy lubricants. This value of H can be identified 
with the structural strength of the bushes saturated with the 
lubricating oil, considered the basic parameter of the durability 
of porous slide bearings. 


POSSIBILITY OF TRANSFER OF THE 
MEASURED RESULTS TO A REAL OBJECT 


The rate of the conversion of the porous slide bearing from 
the state of operating ability to partial ability, or disability, can 
be presented using the Weibull, normal, or another distribution. 
Therefore transferring the experimental data from the physical 
model to a real object is possible, for instance by creating the 
interval of porous slide bearing technical efficiency tolerance 
[8] depending on the bush structure parameter which determines 
the durability of porous slide bearings. This can be done in the 
following steps: 

e measuring the strength of bonds of the greasy lubricant 
having the characteristics correlated (via acoustic emission, 
for instance [7]) with the strength of the porous bearing bush 
structure with respect to the operating time qt, (the values 
will be within the interval u’—u” and their distribution will 
be close to normal — Fig. 9); 


fu) 


T 


Fig. 9. Interval of porous slide bearing technical efficiency tolerance: 
t — time of porous slide bearing operation, u — strength of porous slide bearing 
bush structure, u „— bush structure strength limit, f(u) — probability density 
of the random variable u, f(t) — probability density of the random variable t 


e determining the time period of full operating ability of the 
porous slide bearings until time t, with the tolerance ranging 


Table 1. Evaluated values of diagnostic parameters 


Results of measurements for n = 500 rev/min 


Seizing load P[N] 


Graphite grease 


Seizing load limit P [N] | Strength of bonds H[Nm] 


Grease LT43 


Composition of greases 


Results of measurements for n = 1440 rev/min 


Seizing load P[N] 


Graphite grease 


Seizing load limit P [N] | Strength of bonds H[Nm] 


Grease LT43 


Composition of greases 
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within T’ - t” for the known limit of the strength of bonds 
of the greasy lubricant (the value taken from experimental 
investigations) - Fig. 9. 


FINAL REMARKS AND CONCLUSIONS 


- Although the physical model of a porous slide bearing, 
presented in the article as the friction pair of the apparatus 
T-02, cannot be treated as fully adequate to real porous slide 
bearings, the obtained results of the tests indicate that the 
application of the here presented method for examining the 
strength of the porous bush structure is fully justified. 

- The presented theoretical and experimental analysis of 
the physical model shows that the strength of bonds of 
the greasy lubricants, i.e. their lubricating ability can be 
considered one of parameters of porous bearing bush 
structure. The analysis of lubricating ability of the greasy 
lubricants determines the area of applications of the porous 
slide bearings with respect to their lubricating ability, i.e. 
durability. 
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Influence of local bush wear on properties 
of water lubricated marine stern tube bearings 


Wojciech Litwin, Ph. D. 
Gdansk University of Technology 


ABSTRACT 


In recent years more and more frequently can be met ships in which propeller shaft 
water-lubricated polymer bearings have been applied. It results from their simplicity 
and associated relatively low initial and operational cost as compared with the complex 
classical sealed system based on oil-lubricated bearings. It is worth mentioning that the 
water-lubricated bearings are also environmentally friendly as no risk of pollution with 
lubricating oil used in classical systems, is involved. Design procedure of bearings in 
question based on materials available from their producers is relatively simple. However 


the calculating of such bearings creates many problems. For calculating the hydrodynamic lubricating 
theory is as a rule used and elastic bush flexibility taken into account. In the case of greater loads when 
lubricating film is very thin to achieve a proper solution is very difficult. 
For this reason is proposed a novel method in which bush local wear effect can be accounted for. Such 
phenomenon has been observed during experimental tests of heavy-loaded water-lubricated slide bearings. 
With the use of the applied calculation algorithm bush local wear close to its edges can be taken into 
account, and - in consequence - a greater bearing capacity can be achieved. 


Keywords: bearing systems; water-lubricated bearing; propeller shaft bearing systems 


INTRODUCTION 


In recent years water-lubricated bearings fitted with 
polymer bushes have been more and more commonly applied 
to shipbuilding, hydro power industry and water pumps. 
Such solution has many merits among which lack of risk of 
contamination of the environment by lubricants is crucial one. 

Recently more and more sliding polymers capable of 
operating in water have been appeared on the market. Their 
price has been dropping that makes such solution more and more 
attractive. Moreover, many data on successive applications of 
water - lubricated bearings have appeared. Their producers have 
used the issues for advertising their products whose life time is 
estimated to about thirty years. However as results from practice 
such bearings behave in very different ways during operation 
and in the past many dangerous breakdowns resulting in ship 
docking or premature repairing of water turbine, happened. 

The bearing lubricated with water, a low viscosity liquid, 
shows many limitations. Doubtlessly, it develops always 
a smaller hydrodynamic capacity as compared with oil- 
lubricated bearing of the same size. In addition, it shows form 
instability of polymer bush which, due to shrinkage mounting 
and water soaking, is capable of suffering large deformations. 
There are known many cases of propulsion shaft or rudderstock 
blocking due to polymer bush swelling. In order to prevent 
shaft seizing bush producers often recommend to increase 
bearing clearance, that results in lowering their hydrodynamic 
features. A less dangerous consequence of bush expansion 
during shrinkage mounting and subsequent water absorption 
in initial phase of its operation is rather large form deformation 
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of its sliding surface. The effect is especially distinct if lead- 
over longitudinal grooves necessary for heat transferring by 
water from friction zone, are made in the bush. In such case 
significant deformations occur in the vicinity of the grooves as 
bush material is there of the smallest thickness [1]. 

In calculating a water-lubricated bearing an isothermal 
model can be used, that highly facilitates to obtain calculation 
results [2]. In order to achieve reliable results for bearing 
hydrodynamic capacity elastic polymer bush flexibility should 
be accounted for. However it makes a trouble. It turns out that 
to precisely determine shape elasticity value is rather difficult. 
Usually only its approximate value is given. Hence especially 
in the case of more flexible materials calculation results may 
contain a great error, that was confirmed experimentally [3]. 


RESEARCH PROBLEM 


The research problem in question results from the need of 
answering the question as to which extent water-lubricated 
polymer hydrodynamic bearing can be loaded. 

In the Faculty of Ocean Engineering and Ship Technology, 
Gdansk University of Technology, have been conducted for 
several years experimental and theoretical projects aimed at 
solving a set of problems associated with water-lubricated 
bearings. One of the crucial research problems was to determine 
ultimate hydrodynamic load-carrying capacity of bearings of 
the kind [4, 5]. Results of the performed experimental tests have 
been very promising. Polymer bearings of many kinds were 
tested. The shaft of the test stand was of 100 mm diameter, 
and the tested bushes were from 100 to 400 mm in length. The 


shaft rotational speed was varied in the range of 0 + 11 rev/s, 
that made it possible to obtain linear slide velocities below 3.5 
m/s. From the performed tests it results that the bearings have 
operated in the fluid friction range even under the pressure of 
the order of 0.5 MPa (Fig. 1 and 2). 


0.3 
load 0.3 MPa 
load 0.4 MPa 
0.1 load 0.5 MPa 
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Fig. 1. Diagram of the measured friction torque converted to value 
of coefficient of friction, under the load up to 0.5MPa and 
at shaft rotational speed in the range of 1 + 11 rev/s 
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Fig. 2. The pressure distribution in bearing lubricating film, measured on 
the test stand at the shaft rotational speed of 11 rev/s and under the load 
range from 0.3 to 0.5 MPa: a) in circumferential direction within the range 
of 0 = 22, b) in longitudinal direction within the range of 0 + 300 mm 


A serious problem is met in conducting calculations for 
the similar bearing by applying the hydrodynamic lubrication 
theory, (EHL), to bush elastic material. Namely, below certain 
minimum film thickness which, in the case in question, amounts 
to about 5 to 7 um for elastic materials, it was not possible to 
obtain a convergent solution. The iteration between a module 
which calculates pressure distribution in lubricating gap and 
a module which calculates elastic deformation field of bush 
surface has not led to any convergent solution (Fig. 3) 
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Fig. 3. An attempt to calculating hydrodynamic capacity of a heavy-loaded 
bearing fitted with elastic bush, a, b) for film thickness values of over 5um 
and c) for a smaller value - in this case the iteration diverges 
and no solution is available. 


Attempts to solving the problem which consists in 
introduction of the so-called relaxation factor, did not give 
any satisfying result. In practice in the case in question it was 
not managed to achieve a stable solution for lubricating films 
thinner than 5 + 7um. The calculations performed for the 
minimum acceptable lubricating gap of Sum showed a much 
lower value of bearing capacity as compared with results of 
experimental tests (Fig. 4). 
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Fig. 4. The calculated maximum potential hydrodynamic capacity 
of water-lubricated bearing in function of bush material, 
bearing diameter clearance; shaft speed 11 rev/s 


However it is possible to reach certain conclusions by 
analyzing results of the calculated pressure distributions and 
bush elastic deformations (Fig. 5 and 6). As turns out, elastic 
deformation value may exceed that assumed for calculations of 
minimum film thickness amounting to Sum. In practice an elastic 
sag of bush results in dropping value of hydrodynamic pressure 
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in lubricating film, that detrimentally affects bearing capacity. 
As appears, the zones close to bush edges are not elastically 
deformed due to hydrodynamic pressure action. Analyzing the 
experimentally tested bushes one observed significant wear in 
the zones, which is visible because sliding surface fast became 
polished in the places. 
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Fig. 5. The calculated pressure distribution in lubricating film of water- 
lubricated polymer bearing, (P1), in function of bush length 
(of 1 +200 mm) and bush developed circumference (of 0 + 27 rad) 
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Fig. 6. The calculated distribution of elastic deformation of polymer bush 
of water-lubricated bearing, (P1), in function of bush length 
(of 1 + 200 mm) and bush developed circumference (of 0 + 27 rad) 


PRESENTATION OF THE PROPOSED 
CALCULATION METHOD 


During searching subject-matter literature sources dealing 
with slide bearing systems, has been found several papers 
whose authors undertook the task to predict growing wear [6 
+ 10]. Unfortunately the papers did not help in solving the 
wear problem of water-lubricated bearing and its effect on 
bearing capacity. 

Basing on the earlier conducted projects one decided to 
perform calculations for the case when local wear occurs close 
to bush edges. To make it possible the classical algorithm of 
EHL method has been modified. 

It was assumed that the minimum lubricating gap height was 
equal to Sum. If due to application of greater and greater load 
the lubricating gap decreases and its height becomes smaller 
than the assumed minimum then wear will occur. As a result 
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of the process bush geometry undergoes permanent change. 
The schematic diagram of the relevant calculation software is 
shown in Fig. 7. 

It was assumed that the method in question has to serve 
for calculating slide bearing capacity with accounting for edge 
wear. However if the procedure of application of greater and 
greater load is not terminated then the model will approach to 
the ridging effect reported in literature sources, when the journal 
rubs out a longitudinal groove in the bush. In the author’s 
opinion the effect will not occur in the case of ship shaft line 
bearing as the shaft operates with various rotational speeds and 
the bearing’s loading results from many different factors, e. g. 
ship operation in heavy weather conditions. Therefore working 
position of the shaft against the bush changes in operation and 
bush ovalization appears after long period of service (Fig. 8). 


Fig. 8. Heavy-used water-lubricated slide bearing in which distinct bush 
ovalization took place and consequently the shaft altered its position down 
by the value Ow; ødz external bush diameter, adc initial bush internal 
diameter, adz worn area diameter 


THE TESTED BEARING AND ITS 
SIMULATED WORKING CONDITIONS 


The calculations were performed for two kinds of bush 
material of different values of module of shape elasticity. The 
analysis was performed for the polymer of shape elasticity 
module equal to about 2500 MPa and the rigid composite of 
shape elasticity module equal to about 4500 MPa (Tab. 1). As 
the lubricating medium fresh water of 5°C temperature and 
viscosity of about 1.55 mPa-s, was assumed. The calculations 
were conducted for the shaft rotational speed of 11 rev/s. 


RESULTS OF THE CALCULATIONS 


Results of the calculations of hydrodynamic pressure 
distributions in lubricating film, and bush elastic deformations 
resulting from interaction between pressure field and wear field, 
are presented in Fig. 9 + 12. The developed surface diagrams 
were made for the lower part of the bush of the analyzed bearing 
(for circumference arc changeable from 0.52 to 1.5 m) and the 
bush length from 0 to 200 mm. In order to make interpretation 
of the diagrams easier, values of elastic deformation and 
wear are shown within the same range, i. e. from 0 to 

20 um. Values of ultimate relative journal-bush eccentricity 
as well as results of the calculated bearing capacity are shown 
in Tab. 1. 


START 
Data input, preliminary calculations 
Bush size, shaft speed, water viscosity, 
shaft eccentricity ratio Esan, Emax maximum expected load Lmax, Minimum accepted film thickness hmin, 


Calculations of bearing interspace geometry 


Calculations of pressure field 
Procedure based on hydrodynamic theory of 
lubrication — Reynolds equation 
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Fig. 7. The schematic diagram of the proposed calculation method 
of elastic hydrodynamic lubrication with accounted for bush local wear, (marked EHL+W) 
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Fig. 9. The distribution of the hydrodynamic pressure in lubricating film for Fig. 10. The distribution of elastic deformation for the bush of the elasticity 
the bearing fitted with elastic bush of the elasticity module of 2500 MPa, module of 2500 MPa, (P2), at the shaft rotational speed of 11 rev/s 


(P2), at the shaft rotational speed of 11 rev/s 
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Tab. 1. Input data and results of the calculations for two kinds of the analyzed polymer bush of different values of shape elasticity module 


Bush 
material 


type 


Module of 
elasticity 


Bearing code 


bush wear [mm] 


circumferential direction [rad] 1.50 
Fig. 11. The wear field for the bush of the elasticity module of 2500 
MPa,(P2), at the shaft rotational speed of 11 rev/s; the wear occurred close 
to the bearings edges only 


— a 
— =y 


bush wear [mm] 


circumferential direction [rad] 1 5r 


Fig. 12. The wear field for the bush of the elasticity module of 2500 MPa, 
(P2),; the wear occurred also along the whole breadth of the bearing; this is 
the case when the calculations have been discontinued; the shaft rotational 
speed of 11 rev/s; the calculated bearing capacity of 6500 N 


CONCLUSIONS 


The proposed calculation method based on the EHL model 
extended by a module which accounts for bush wear, made it 
possible to obtain calculation results for the conditions in which 
the classical model ceases to be effective. 

In the author’s opinion the proposed calculation method is 
reliable as the phenomenon of local seizing close to bush edges 
was observed in the experimentally tested bearings. 

The performed calculations demonstrated that the accounting 
for bush local wear makes it possible to increase theoretical 
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Bearing 
clearance/shaft 
diameter/bush 
thickness [mm] 


0.3/100/12 


Shaft 
eccentricity 
ratio (£) 


Calculations 
results — maximum 
hydrodynamic capacity 


Maximum 
wear [um] 


bearing capacity by 7% in the case of the more rigid composite 
material and by nearly 15% for the less rigid one (polymer). 

The calculated bearing hydrodynamic capacity with 
accounted for bush local wear reaches the value of 0.316 MPa. 
On the basis of experimental tests it is possible to conclude that 
the maximum hydrodynamic capacity of bearings of the kind 
may reach even the value of 0.5 MPa. Therefore it should be 
considered why not to carry out calculations for smaller values 
of lubricating film thickness and greater values of bush wear 
located not only close to its edges but also along the whole 
contact line between shaft and slide bush (Fig. 12). 

The proposed calculation method may be applied also in the 
case of the skewed shaft against bush, i.e. in which the method 
may be especially helpful. 
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Model of thermal comfort 
in the hyperbaric facility 
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ABSTRACT 


The paper discusses the mathematical model of thermal comfort in the hyperbaric facility. Based on human 
thermal balance and thermal comfort conditions the comfort equation for the hyperbaric environment was 
derived. The comfort equation enables to calculate all those combinations of the diver 5 activity, clothing 
(thermal insulation, moisture permeation factor) and environmental variables (temperature, pressure, 
composition, relative humidity and velocity of the breathing gas and mean radiant temperature), which 
would create thermal comfort in the hyperbaric environment. The paper presents also the solution of 
thermal comfort equation for dives up to the depth of g = 590 m, with helium-oxygen and hydrogen-oxygen 
breathing mixture, as well as the comfort diagrams. 


Keywords: comfort temperature; human balance; hyperbaric environment 


INTRODUCTION 


Saturation diving is the diving operation that enables man 
to live and work in the sea-going deep dive systems or to be 
trained in the land-based hyperbaric complexes for a long time. 
Advanced diving technique requires an optimal breathing gas 
to be selected. The use of air as the breathing gas is limited 
because of oxygen toxicity and narcotic effect of nitrogen 
[5, 14, 15]. Breathing mixtures used in deep diving consist 
of oxygen that meets demands for metabolic needs and inert 
gases (helium, nitrogen, hydrogen). Moreover, breathing gas 
contains contaminants generated by the divers during exposure 
and some gaseous additives improving voice communication 
(sulphur hexafluoride, carbon tetrafluoride) or freon used as 
fire extinguishing agent [5]. 

Saturation diving is based on the principle that the tissues 
and blood of a diver’s body absorb inert gases and reach 
saturation point where the pressure of the dissolved gas in the 
blood and tissues is the same as that of the gas in the lungs. The 
diving technique is based on the fact that once the saturation 
point has been reached the diver can safely remain saturated 
for a long period, and the diver’s obligation for decompression 
does not increase with further exposure. It makes deep diving 
operations more efficient. 

Most of the deep diving experience is based on the breathing 
gas containing helium, however the last years have shown that 
hydrogen also can be used [4] as the component of breathing 
mixtures used in saturation diving There were also some trials 


with experimental breathing mixtures containing the other inert 
gases: neon, argon [5]. 

Helium and hydrogen causes certain thermal problems 
related to dangerous cooling of the divers. This is due to the 
high thermal conductivity and specific heat of helium and 
hydrogen. Therefore, comfort temperature in the hyperbaric 
environment increases with the total pressure of gas and should 
be maintained at a level higher than that in atmospheric air 
[8, 9, 11, 14, 15, 16, 17, 18, 19, 20]. 

A see-going deep dive system presented in Fig.l consists 
of the support vessel, deck decompression chamber and diving 
bell. A land-based hyperbaric facility presented in Fig. 2 
includes one or a series of the hyperbaric chambers where 
the divers live and work at the pressure which simulates the 
ambient one at their work site. Dry or wet chambers of the 
land-based hyperbaric facility are connected by locks. The 
chambers can be maintained at the same pressure or different 
one depending on the goals of hyperbaric tests (medical tests, 
equipment tests, etc). 

One of the major tasks of the life support systems of 
the hyperbaric facilities is to create the environment which 
maintains the safe level of oxygen partial pressure, required 
level of contaminants and the thermal comfort for the divers 
[15]. The thermal comfort is defined by ASHRAE standard 
55-66 [2], as a state of mind that express satisfaction with the 
surrounding environment. 

There are many variables which influence the thermal 
comfort condition in hyperbaric environment: pressure, 
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Fig. 1. Sea-going deep dive hyperbaric system: 1) support vessel, 
2) deck decompression chamber, 3) diving bell, 4) diver 


> 
3 


Fig. 2. Land-based hyperbaric complex: 1) dry hyperbaric chamber, 
2) lock, 3) ,,wet-pot” hyperbaric chamber 


composition, temperature, humidity and relative velocity of 
the breathing mixture, mean radiant temperature, activity level 
and thermal insulation of the diver’s clothing. It is impossible 


1. DIVING OPERATION 


diving facility 

depth of diving 
number of the divers 
duration of exposition 
decompression 


fatigue 
acclimatisation 
health 
type of diet 
noise 
failure of equipment 


to consider the effect of any of the above mentioned variables 
independently. Thermal comfort can be reached by many 
different combinations of the environmental variables, however 
there are also some irrational (human) factors (stress, fatigue, 
etc.) which should be taken into consideration. Fig. 3 presents 
the factors (grouped in six blocks - pieces of thermal comfort 
puzzle supplemented with block of irrational factors.) which 
affect thermal comfort in the hyperbaric facility. 

Not many earlier comfort studies for the hyperbaric 
environment have taken all the above mentioned variables 
into account. However, few studies have been performed on 
men to ascertain the comfort conditions experimentally and 
theoretically [3, 8, 9, 11, 13, 16, 18, 19, 20]. 

In order to determine the comfort temperature in the 
hyperbaric environment the comfort equation based on the 
mathematical model of human body thermal exchange, was 
developed. Mathematical modelling of thermal comfort enables 
to assure the comfort temperature by taking into account the 
most important variables which affect thermal comfort in the 
hyperbaric facility: 


t =f(P, X; C sAN Ps Åmo Leas tr Ps V) (1) 


THERMAL COMFORT EQUATION FOR 
THE HYPERBARIC FACILITY 


The overall steady- state thermal balance of the human body 
can be expressed [2] with the following equation: 

The equation (2) expresses that the flux of metabolic heat 
production, Q,,, minus the flux of heat loss by evaporation 
from the skin, Q., and that by respiration, Q., is equal to 
the flux of heat conducted through the clothing, Q , and 
that dissipated from the outer surface of the clothing by 


convection, Q„ and radiation, Q.. Figure 3 shows the overall 
thermal balance of the human body. As seen from Fig. 3 the flux 


Fig. 3. Factors which influence thermal comfort in the hyperbaric facility. 
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Qm 


Qj 


Qu 


Fig. 4. Overall thermal balance of the human body 


ofheat loss by respiration, Q » and by evaporation from the skin, 
Q,, are not driven by the temperature difference directly. 

The flux of metabolic heat released in oxidation processes in 
human body is partly converted into the external power N and 
the internal body heat flux Q, hence thermo-regulatory system 
maintains the core temperature t, approximately constant: 


Qn =m Âp =Q+N (3) 
External mechanical efficiency is defined as [2]: 
N 
NEK (4) 


Qui 


By introducing the equation (4) to equation (3) the following 
is obtained: 


Q=Q,,0-n4a) (5) 

The flux of heat loss by respiration consists of the flux of 
dry respiration heat loss, Qop and the latent respiration heat 
loss Q: 


ou’ 
Q, = Qij + Qou (6) 
The flux of heat loss by evaporation consists of flux of 


heat loss by water vapour diffusion through skin and sweat 
evaporation from the skin: 


Q. = Qa + Qua (7) 

The heat loss by sweat evaporation is affected by skin 
temperature and the activity level. The skin assists in thermal 
homeostasis. The blood vessels and sweat glands located in 
the skin are activated due to thermal environmental conditions. 
Vasodilatation which allows increased blood to flow through the 
artery, sweat secretion and evaporation are the most important 
ways by which man attempts to lose excess of body heat. In 
cold conditions sweat is not secreted, muscles under the surface 
of the skin contract lifting the hair follicle upright and cause 
“goose bumps”. Vasoconstriction enables to direct blood away 
from the skin and towards the warmer core of the body. This 
prevents blood against losing heat to the surroundings and also 
prevents the core temperature against decreasing. 

The transfer of heat from skin to the outer surface of clothed 
body is rather complex as it consists of internal convection and 
radiation in material pores and conduction through material 
of the clothing itself, [2]. Flux of heat transferred from the 
skin to the outer surface of clothed body depends upon the 
total heat resistance to transfer from skin to the outer surface, 
skin temperature and temperature of the outer surface of the 
clothing: 


Q = t, = ty) i Apu 
R 
The total heat resistance to transfer from skin to the outer 
surface of the clothed body is given as follows [2]: 
0.155 7 
cla 
Ry 
The flux of heat dissipated by convection from the outer 
surface of the clothing, Q o and radiation, 0; depend upon 
ambient temperature and mean radiant temperature, emissivity 
of skin and clothing surface and thermal properties of the 
environment. Satisfaction of the heat balance equation (2) is the 
first sufficient condition for the diver’s thermal comfort under 
long exposure to the hyperbaric environment. At a given activity 
level the skin temperature and the heat loss by evaporation are 
deemed the only physiological variables which influence the 
heat balance equation (2). Sensation of thermal comfort has 
been related to magnitude of the two variables: 


t =f), a<t, <b 


(8) 


cl 


(10) 
de= fy (m) ? c<, <d (11) 


Man is in thermal comfort when the mean skin temperature 
and sweat secretion is maintained within the narrow limits 
(a, b) and (c, d), respectively, determined by the physiologists. 
According to the mathematical comfort model, [2], the 
equations (10) and (11) represent the second and third basic 
conditions for thermal comfort. The skin temperature for 
persons in thermal comfort, experimentally investigated in 
atmospheric air, can be expressed [2] as follows: 


t,=35.7 -0.026-q,, -(l— ny) (12) 


Under the assumption that the skin temperature in hyperbaric 
environment is nearly the same [3] as in atmospheric air it is 
more convenient to express the second comfort condition (10) 
in function of the oxygen consumption, external mechanical 
efficiency and type of diet: 


t, = 35.7 -0.032-C-Vo, :(1-n4) (13) 


Average value of the coefficient C which describes the 
utilization of 1 cm? volume of oxygen can vary from C = 18.85 
to 20.95 J/cm? depending on the proportion of carbohydrate, 
protein and fat being oxidized in the diver’s organism [3]. The 
second basic comfort condition [2] enables to determine heat 
loss by evaporation of sweat secretion. It can be expressed for 
the hyperbaric environment as follows: 
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Q. 
A 


Area of the nude body surface, i.e. DuBois area, is expressed 
by the following equation [2]: 


A. =0.203 - m?425 « 725 
Du . 


Skin wetness is defined [1] as follows: 


d. = =w- p-r: Eas- P) (14) 


Du 


(15) 


w=- Q. = 0.02+0.4[l- exp(— 0.01034 -d „+ 0.6)] (16) 
emax 
The clothing moisture permeation factor F, describes the 
„Cooling efficiency” of sweating on the surface for a clothed 
human body [10]. The permeation efficiency factor F „has been 
derived theoretically for the different kinds of the hyperbaric 
breathing mixtures [8]: 


=] 
pla af} OSSR TAB] Gy 


dm w-m 


By substituting all the terms of heat loss derived in [8], and 
thermal comfort basic conditions into the heat balance equation 
(2) the following is obtained: 


4.58 +3.46-v 
—0.8895-10" -åm T- (P, -PP)+ 


Åm (a =a) 


—4.1-10° -g,,-¢, P R” (ta —t)+ 


=w- Per Fa’ s -O Pe) = (18) 
à- (t ta) 
4.03 -10 > -Ioa a ‘(tat 


+4-10° -f, (T4 -T’) 


The comfort equation derived for the hyperbaric environment 
can be expressed in function of breathing mixture pressure, 
composition and thermal properties (specific heat, thermal 
conductivity, viscosity and density), diver’s activity, thermal 
resistance of the clothing, gas temperature, mean radiant 
temperature, relative humidity and velocity of the breathing 
mixture, (18). 

The comfort equation (18) makes it possible to calculate 
all combinations of the above specified variables which are 
responsible for thermal comfort in the hyperbaric chamber. 
Moreover, it is possible to determine the comfort temperature 
inside the hyperbaric facility and all terms of the heat dissipated 
from the human body exposed to the hyperbaric environment. 
It’s worth to accentuate that thermal comfort is felt very 
individually. 


THE RESULTS 


The comfort equation was solved under the following 

assumptions: 

e the depth of simulated diving: g = 90, 190, 290, 390, 490, 
590 m; 

e the breathing gas: oxygen + helium, (HELIOX), oxygen 
+hydrogen (HYDROX); 

e the relative humidity of the breathing gas: ọ = 0.4, 0.5, 0.6, 
0.7, 0.8, 0.9, 1 at each depth; 

e the relative velocity of the breathing gas: v= 0.1, 0.15, 0.2, 
0.25, 0.3 m/s; 
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e the mean radiant temperature: t,=0, 10, 20, 30, 40, 50 °C; 

e the Du Bois area: A, =1.8 m?; 

e the diver was clothed in a clothing of the total heat transfer 
resistance from skin to the outer surface of 1, = 0.2 + 1.1 
clo, clothing area factor f, = 0, 8, 1 clo = 0.155 m’K/W; 

e the activity levels: the sedendary activity q,, = 58 W/m’, 
medium work activity q,, = 90 W/m’, intensive work 
activity q,, = 120 W/m’. 


The comfort equation has been solved for relevant 
combinations of the above specified variables and the results 
have been presented in the form of the diagrams. The curves 
in Fig. 5 and Fig. 6 (gas pressure versus comfort temperature 
with gas velocity as a parameter) represent the comfort lines i.e. 
lines drawn through points which satisfy the equation (17) and 
thus provide the thermal comfort in HELIOX and HYDROX 
breathing mixtures under the following assumptions: the 
relative humidity of breathing gas ọ = 0.7, sedendary activity 
level įm = 58 W/m’, clothing thermal resistance I „= 0.37 clo, 
mean radiant temperature 1, = 0.37. 


HELIOX 
33.5 


Comfort temperature, t[°C] 


0.5 1.5 2.5 3.5 4.5 3.5 
Pressure, p [MPa] 


Fig. 5. Comfort lines for HELIOX (pressure versus comfort temperature 
with parameter of relative gas velocity) under the following assumptions: q, 
= 58 Wim’, I „= 0.37 clo, p = 0.4, t= 30 °C 
33.6 
33.2 
32.8 


32.4 


32.0 


Comfort temperature, t|°C] 


31.6 
0.5 1.5 2.5 3.5 4.5 5.5 
Pressure, p [MPa] 
Fig. 6. Comfort lines for HYDROX (pressure versus comfort temperature 
with relative gas velocity as a parameter) under the following assumptions: 
Am = 58 Win, I, = 0.37 clo, o = 0.4, t, = 30 °C 


From Fig. 5 and Fig. 6 it is evident that the total pressure 
of the breathing gas exerts a large effect on thermal comfort. 
To compensate body heat losses and achieve thermal comfort 
within the gas pressure range of p = 0.98 + 5.88 MPa and gas 
relative velocity range of v = 0 + 0.3 m/s it is necessary to 
increase the comfort temperature by At = 2.2 K. As it follows 
from calculations despite pressure range, at p = const the 
increase of relative gas velocity from v = 0 to v= 0.3 m/s should 
be compensated for by the comfort temperature increase of 
At=0.7 K. 

Fig. 7 and Fig. 8 show comfort lines which illustrate the 
relationship between comfort temperature in HELIOX and 
HYDROX atmospheres, mean radiant temperature and relative 


gas velocity under the following assumptions: d „= 58 W/m’, 
La = 0.37 clo, @ = 0.7, t, = 30 °C, p = 0.98 MPa. As it follows 
from Fig. 7 and Fig. 8 at constant relative gas velocity 
increase of mean radiant temperature should be compensated 
by comfort temperature decrease. As seen from Fig. 7 and 
Fig. 8, the comfort lines cross each other at the mean radiant 
temperature t, = 20 °C and comfort temperature t = 34.1 °C 
for HELIOX, and t, = 20 °C and t = 34°C for HYDROX, 
respectively. 


HELIOX 


Comfort temperature, t[°C] 


0 10 20 30 40 50 
Mean radiant temperature, t, [°C] 
Fig. 7. Comfort lines for HELIOX (mean radiant temperature versus 


comfort temperature with relative gas velocity as a parameter) under the 
following assumptions: q,, = 58 W/m, 1, = 0.37 clo, p = 0.7, p = 0.98 MPa 


HYDROX 


Comfort temperature, t[°C] 


0 10 20 30 40 50 
Mean radiant temperature. t, [°C] 
Fig. 8. Comfort lines for HYDROX (mean radiant temperature versus 


comfort temperature with gas velocity as parameter) at the following 
assumptions: In = 58 W/m’, LS 0.37 clo, o = 0.7, p = 0.98 MPa 


The comfort lines cross each other where the temperature 
of the breathing mixture is equal to the mean temperature of 
the clothed body since the convective heat transfer here will 
be zero, independent of the relative velocity. To the right of 
the crossing point the temperature of the outer surface of the 
clothing is higher than the breathing mixture temperature and 
an increase in the relative velocity will therefore require an 
increase of the breathing mixture temperature or changing mean 
radiant temperature in order to maintain thermal comfort.To the 
left of the crossing point the temperature of the outer surface of 
the clothing is lower than the breathing mixture temperature and 
an increase in gas velocity will therefore require a decrease in 
the mixture temperature or changing mean radiant temperature 
to maintain thermal comfort. 

Fig. 9 and Fig. 10 show comfort lines for HELIOX and 
HYDROXxX, which illustrate the relationship between comfort 
temperature, pressure and the total heat transfer resistance 
of clothing under the following assumptions: q,, = 58 W/m’, 
ọ = 0.4, t, = 25 °C, v = 0 m/s. 

From Fig. 9 and Fig. 10 it is evident that for all types of 
the total heat transfer resistance of clothing increase of gas 
pressure should be compensated for by increase of comfort 
temperature. As it follows from Fig. 9 and Fig. 10 at I „= const, 


increase of pressure from p = 0.98 MPa to p = 5.88 MPa should 
be compensated for by increase of comfort temperature by 
At = 0.7 K for HELIOX and by At = 0.6 K for HYDROX. To 
maintain thermal comfort in HELIOX atmosphere of pressure 
p = 2.94 MPa and in HYDROX atmosphere of pressure 
p = 1.96 MPa increase of the total heat transfer resistance 
of clothing from I, = 0.3 clo (light clothing ensemble) to 
ILa = 1.1 clo (heavy working ensemble) should be compensated 
for by decrease of comfort temperature of At = 1.3 K for 
HELIOX and of At = 0.8 K for HYDROX. 


HELIOX 


Comfort temperature, t[°C] 


0.5 1.5 2.5 3:5 4.5 5.5 


Pressure, p [MPa] 


Fig. 9. Comfort lines for HELIOX (pressure versus comfort temperature 
with the total heat transfer resistance of clothing as a parameter) under the 
following assumptions: q,, = 58 Wim’, 9 = 0.4, t, = 25 °C, v = 0 m/s 
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Fig. 10. Comfort lines for HYDROX (pressure versus comfort temperature 
with the total heat transfer resistance of clothing as a parameter) under the 
following assumptions: q,, = 58 Wim’, 9 = 0.4, t, = 25 °C, v = 0 m/s 


Fig. 11 and Fig. 12 show the relationship between comfort 
temperature, pressure and the diver’s activity in HELIOX and 
HYDROX atmosphere under the following assumptions: @ = 
0.6, t, = 30 °C, v= 0 m/s, I „= 0.2 clo. It is evident from Fig. 
11 and Fig. 12 that at the constant activity of the diver, comfort 
temperature increases due to increase of pressure. To maintain 
thermal comfort at p = const and increasing activity of the diver 
it is necessary to reduce comfort temperature. 

HELIOX 
34.0 
33.0 
32.0 + 
31.0 
30.0 
29.0 
28.0 


27.0 


120 W m? 


Comfort temperature, t[°C] 


26.0 
0.5 1.5 2.5 3.5 4.5 $5 
Pressure, p [MPa] 

Fig. 11. Comfort lines for HELIOX (pressure versus comfort temperature 


with the diver activity as a parameter) under the following assumptions: 
L7 0.2 clo, = 0.6, t, = 30 °C, v = 0 m/s 
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HYDROX 


Comfort temperature, t[°C] 
N w 


— o 58 Wim? == q = 90 Wim qn = 120 W m 
0.5 1.5 2.5 3.5 4.5 5.5 
Pressure, p [MPa] 
Fig. 12. Comfort lines for HYDROX (pressure versus comfort temperature 
with the diver activity as a parameter) under the following assumptions: 
la” 0.2 clo, 9 = 0.6, t =30 °C, v = 0 m/s 


As it follows from Fig. 11 and Fig. 12 to maintain thermal 
comfort at the pressure p = 0.981 MPa, change in the diver’s 
activity from q,, = 58 W/m? (sedendary) to q,, = 120 W/m? 
(intensive work) should be compensated for by the comfort 
temperature decrease At = 4.9 K in the case of HELIOX, 
and At = 3.0 K in the case of HYDROX. At the diver’s 
sedendary activity q = 58 W/m? the pressure increase from 
p = 0.98 MPa to p = 5.88 MPa should be compensated for by 
the comfort temperature increase At = 1.5 K both for HELIOX 
and HYDROX. 

Fig. 13 and Fig. 14 show comfort lines for HELIOX and 
HYDROX under the following assumptions: I „„ = 0.37 clo, 
ọ = 0.6, t, = 30 °C, p = 1.96 MPa. As it follows from Fig. 13 
and Fig. 14 at constant activity an increase of relative gas 
velocity should be compensated for by an increase of comfort 
temperature. At v = const an increase of diver’s activity should 
be compensated for by a decrease of comfort temperature. 

HELIOX 
33.0 
32.0 
31.0 
30.0 
29.0 
28.0 


27.0 


Comfort temperature, t[°C] 


0 0.05 0.1 0.15 
Gas relative velocity, v [m/s] 
Fig. 13. Comfort lines for HELIOX (gas velocity versus comfort 
temperature with the diver activity as a parameter) under the following 
assumptions: 1 „= 0.37 clo, o = 0.6, Ls 30 °C, p = 1.96 MPa 
HYDROX 
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0 0.05 0.1 0.15 0.2 0.25 0.3 
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Fig. 14. Comfort lines for HYDROX (relative gas velocity versus comfort 


temperature with the diver activity as a parameter) under the following 
assumptions: 1, = 0.37 clo, o = 0.6, ic 30 °C, p = 1.96 MPa 


As seen from Fig.13 and Fig. 14, to maintain the sedendary 
diver (q,, = 58 W/m’) and working diver (q,, = 120 W/m’) in 
thermal comfort, when relative gas velocity changes from v = 0 
to v=0.3 m/s, it is necessary to increase comfort temperature 
by At = 0.7 K for HELIOX and At = 3.5 K for HYDROX. As 
shown in Fig.13 and Fig. 14 the change in the diver’s activity 
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from q,, = 58 W/m? to dm = 120 W/m”, at constant relative 
gas velocity of v = 0.1 m/s, should be compensated for by the 
temperature decrease At = 4.4 K for HELIOX and At = 3.5 K 
for HYDROX. 

Fig. 15 shows the relationship between comfort temperature, 
pressure and relative humidity of HY DROX under the following 
assumptions: q,, = 58 W/m’, I, = 0.37clo, v = 0 m/s, mean 
radiant temperature t, = 25 °C. 

As it follows from Fig. 15 at constant gas relative humidity 
an increase in pressure will increase comfort temperature. At 
the pressure p = 0.98 MPa a change in HYDROX relative 
humidity from ọ = 0.4 to ọ = 1 should be compensated for by 
the comfort temperature decrease At = 0.6 K. 

At higher pressure of breathing mixture comfort lines 
are convergent and asymptotically approach the comfort 
temperature t = 33.6 °C. It means that there is no significant 
influence of breathing mixture relative humidity on comfort 
temperature. As it follows from the obtained results the 
influence of the breathing gas relative humidity on the 
magnitude of comfort temperature is not significant however 
very important for sweat evaporation and comfort sensation. 
Poor evaporative rate may result in a diver’s feeling of high 
humidity (discomfort) and will prevent from evaporative 
cooling of the skin. 

HYDROX 


Comfort temperature, t[°C] 
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Fig. 15. Comfort lines for Hydrox (gas pressure versus comfort temperature 
with the relative humidity as a parameter) under the following assumptions: 
In = 58 W/m, I „= 0.37clo, v = 0 m/s, t, = 25 °C, p = 0.98 MPa. 


By using the comfort equation it is possible, for any type of 
the breathing mixture (composition, pressure), type of clothing 
and type of the diver’s activity, to calculate all reasonable 
combinations of the temperature, relative humidity and velocity 
of the breathing gas and mean radiant temperature. The water 
vapour transport coefficient in the hyperbaric environment 
can be computed from the experimental relationship for 
vapour transport coefficient in helium-oxygen mixture [4,5] or 
theoretically [8]. The results of the experiments demonstrated 
a great influence of the total pressure on the water vapour 
transport coefficient in the helium rich atmosphere.To make 
an assessment of the thermal comfort equation derived for 
the hyperbaric facility, its theoretical and an experimental 
verification is recommended. Fig. 16 shows comparison 
between the theoretical results obtained from the thermal 
comfort equation and the experimental results published in 
[12, 16, 18]. 

The theoretical solution of the thermal comfort equation, 
shown in Fig. 16, was obtained under the following assumptions: 
dm = 58 W/m’, I, =0.37clo, v = 0 m/s, ọ = 0.4, t, = 30 °C. As 
it follows from Fig. 16 the theoretical results are in compliance 
with those previously published in [16] for HELIOX, and 
are higher approximately by At = 2 K than the temperature 
predicted in [12, 13, 19]. Comparison of the results is 
difficult. 
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Fig. 16. Comparison between the theoretical results obtained 
from the thermal comfort equation and the experimental results 
published in [12, 16, 18] 


After experimental verification, the mathematical model 


of thermal comfort can serve as a very useful mathematical 
tool in designing the environmental life support systems 
of the hyperbaric facilities. The derived thermal comfort 
equation enables to simulate different thermal conditions in the 
hyperbaric facility with taking into account the most important 
variables which influence thermal comfort. 


CONCLUSIONS 


On the basis of the mathematical model of human body 
thermal exchange and two basic thermal comfort conditions, 
the comfort equation for hyperbaric environment was 
derived 

The comfort equation enables to calculate all those 
combinations of environmental variables which would 
create thermal comfort in the hyperbaric facility, 

By using the comfort equation it is possible, for any type 
of the breathing mixture (composition, pressure), type of 
clothing and the diver’s activity, to calculate all reasonable 
combinations of the gas temperature, relative humidity and 
velocity and mean radiant temperature, 

The comfort equation was solved for the hyperbaric 
breathing mixtures: HELIOX and HYDROX in the 
assumed range of the parameters in diving up to the depth of 
g=590 m, 

As follows from the solution of the comfort equation, the 
required comfort temperature in the hyperbaric facilities 
increases progressively along with increasing pressure 
of HELIOX and HYDROX breathing mixtures, and 
the comfort temperatures are similar in values for both 
atmospheres, 

The results obtained under the assumptions: ġ „= 58 W/m’, 
La = 0.37clo, v = 0 m/s, ọ = 0.4, t, = 30 °C for the pressure 
range of p = 0.98 + 5.88 MPa, have shown the comfort 
temperature range of t = 31.6 + 33.6 °C to be nearly the 
same for HELIOX and HYDROX. 

The theoretical results found for HELIOX are comparable 
with those previously experimentally determined and 
reported in [12, 13, 16, 18]. There are no available data 
concerning the relationship between comfort temperature 
and pressure in H,-O, mixture. 

The presented paper has concluded the first stage of the 
research aimed at mathematical modelling of thermal 
comfort in hyperbaric facility, and the second stage should be 
aimed at experimental verification of the derived model, 
After experimental verification, the mathematical model of 
thermal comfort can serve as a very useful mathematical 


tool in designing and operating the environmental life 
support systems of hyperbaric facilities, as it enables to 
simulate different thermal conditions in the hyperbaric 
facility. 


NOMENCLATURE 


Na 


— Du Bois area, human body surface area, [m°], 

— specific heat of the breathing mixture, [J/kg K], 

— diffusion coefficient of water vapour into breathing 
mixture, [m7/s], 

— ratio of the surface area of the clothed body and the 
surface area of the nude body, 

— moisture permeation factor, 

— thermal resistance to heat transfer from skin to outer 
surface of the clothed body, [clo], 

— diffusive resistance for water vapour through the clothing, 
[m’sPa/kg], 

— diffusive resistance for water vapour through the breathing 
mixture, [m’sPa/kg], 

— external power, [N], 

— total pressure of the breathing mixture, [Pa], 

saturated water vapour pressure at skin (surface) 

temperature and ambient breathing gas temperature, [Pa], 

respectively 

— diver’s activity level, metabolic rate, [W/m’], 

— flux of internal heat produced by the body, [W], 

— flux of heat loss resulting from water vapour diffusion 
through the skin, [W], 

— flux of evaporative heat loss, [W], 

— flux of heat loss resulting from sweat evaporation from the 
skin, [W], 

— flux of evaporative heat loss at skin wetness w = 1, 

— flux of metabolic heat released in oxidation processes, 

[W], 

— flux of respiration heat loss, [W], 

— flux of dry respiration heat loss, [W] 

flux of latent respiration heat loss, [W] 

— flux of convective heat dissipated from the outer surface 

of the clothed body, [W], 

— flux of heat loss from the skin to the outer surface of the 

clothed body, [W], 

— flux of heat loss resulting from radiation from the outer 
surface of the clothed body, [W], 

— specific latent heat of vaporisation, [J/kg], 

gas constant of the breathing mixture, water, [J/kg K], 

— comfort temperature, [°C, 


T, — temperature of the outer surface of the clothing, [°C, K], 


— temperature of expired gas, [°C] 

— mean radiant temperature, [°C, K], 

— temperature of the skin surface, [°C, K], 

— skin wetness, 

— relative velocity of the breathing mixture, [m/s], 

— oxygen consumption, [cm*/min], 

— molar fraction of breathing gas component, 

— convective heat transfer coefficient, [W/mK] 

— water vapour transport coefficient in the hyperbaric 
environment, [kg/m’sPa], 

— relative humidity of the breathing mixture, 

thermal conductivity of the breathing mixture and 

atmospheric air, [W/mK], respectively, 

— viscosity of the breathing mixture, [kg/ms] 

— mechanical efficiency of diver. 
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problems of manoeuvrability (including human effect, harbour and waterway design) are successfully realised at 
the Centre. 

The Foundation possesses ISO 9001 quality certificate. 


Why training on ship handling? 


The safe handling of ships depends on many factors - on ship's manoeuvring characteristics, human factor (operator 
experience and skill, his behaviour in stressed situation, etc.), actual environmental conditions, and degree of water 
area restriction. 


Results of analysis of CRG (collisions, rammings and groundings) casualties show that in one third of all the 
human error is involved, and the same amount of CRG casualties is attributed to the poor controllability of ships. 
Training on ship handling is largely recommended by IMO as one of the most effective method for improving the 
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GDANSK UNIVERSITY OF TECHNOLOGY 


is the oldest and largest scientific and technological academic institution in 
the Pomeranian region. The history of Gdansk University of Technology is 
marked by two basic dates, namely: October 6, 1904 and May 24, 1945. 


The first date is connected with the beginning of the technical education at 
academic level in Gdansk. The second date is connected with establishing of 
Gdansk University of Technology, Polish state academic university. Gdansk 
University of Technology employ 2,500 staff, 1,200 whom are academics. 
The number of students approximates 20,000, most of them studying 
full-time. Their career choices vary from Architecture to Business and 
Management, from Mathematics and Computer Science to Biotechnology 
and Environmental Engineering, from Applied Chemistry to Geodesics and 
Transport, from Ocean Engineering to Mechanical Engineering and Ship 
Technology, from Civil Engineering to Telecommunication, Electrical and 
Control Engineering. Their life goals, however, are much the same - to 
meet the challenge of the changing world. The educational opportunities 
offered by our faculties are much wider than those of other Polish Technical 
universities, and the scientific research areas include all of 21st Century 
technology. We are one of the best schools in Poland and one of the best 
known schools in Europe — one that educates specialists excelling in the 
programming technology and computer methods used in solving complicated 

scientific, engineering, organizational and economic problems. 


THE FACULTY OF OCEAN ENGINEERING 
AND SHIP TECHNOLOGY 


The Faculty of Ocean Engineering and Ship Technology (FOEST) as 
the only faculty in Poland since the beginning of 1945 has continuously 
been educating engineers and doctors in the field of Naval Architecture 

and Marine Technology. 


The educational and training activities of FOEST are supported by 
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international organizations and associations, as well as participation in 
scientific conferences and symposia. Hosting young scientists and students 

from different countries is also a usual practice in FOEST. 


The activities of Faculty departments are related to: mechanics and 
strength of structures, hydromechanics, manufacturing, materials and system 
quality, power plants, equipment and systems of automatic control, mostly 

in shipbuilding, marine engineering and energetic systems. 


FOEST is a member of such organizations like WEGEMT; The 
Association of Polish Maritime Industries and the co-operation between 
Nordic Maritime Universities and Det Norske Veritas. The intensive teaching 
is complemented and supported by extensive research activities, the core 
of which is performed in close collaboration between FOEST staff and 
industry. We take great care to ensure that the applied research meet both 
the long term and short term needs of Polish maritime industry. FOEST 
collaborates with almost all Polish shipyards. Close links are maintained 
with other research organizations and research institutions supporting the 
Polish maritime industry, such as Ship Design and Research Centre and 
Polish Register of Shipping, where several members of the Faculty are also 

members of the Technical Board. 


The Faculty of Ocean Engineering and Ship Technology is a unique 
academic structure, which possesses numerous highly qualified and 
experienced staff in all above mentioned specific research areas. Moreover, 
the staff is used to effective co-operation and exchange of ideas between 
specialists of different detailed areas. This enables a more integrated and 
comprehensive treatment of research and practical problems encountered 
in such a complicated field of activity as naval architecture, shipbuilding 

and marine engineering. 


The staff of the Faculty has strong international links worldwide, being 
members or cooperating with international organizations like International 
Maritime Organization IMO, International Towing Tank Conference ITTC, 
International Ship and Offshore Structures Congress ISSC, International 
Conference on Practical Design of Ship and other floating Structures PRADS 

just to name a few. 
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